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1.   Introduction 

One  of  the  significant  contributions  to  controlled 
thermonuclear  research  in  recent  years  has  been  the 
experimental  work  of  the  Russians  in  developing  their 
Tokamak  machine.   A  useful  mathematical  model   of  the 
Tokamak  can  be  obtained  from  magnetohydrodynamics  in  the 
form  of  a  free  boundary  problem  in  potential  theory. 
Different  methods  have  been  developed  to  compute  equili- 
briiom  configurations  for  cases  which  possess  some  type  of 
symmetry,  for  example,  axial  or  helical  symmetry  [3,5,7]. 

More  important  than  the  computation  of  equilibrium 
configurations  is  the  determination  of  their  stability. 
At  the  present  time,  various  devices  have  been  built,  or 
are  in  the  process  of  being  built.   Some  of  them  are  axially 
symmetric  devices  with  different  kinds  of  cross-sections; 
others  are  nonaxially   symmetric,  as  for  example   those 
with  helical  windings. 

This  raises  the  question  of  which  is  the  best 
configuration  from  the  point  of  view  of  magnetohydro- 
dynamic  stability.   In  order  to  answer  this  question,  it  is 
important  that  we  formulate  a  model  in  which  we  allow  for 
three  independent  space  coordinates . 

The  subject  of  this  paper  is  the  formulation  and 
implementation  of  a  numerical  method  to  compute  configura- 
tions of  magnetohydrodynamic  equilibrium  in  three- 
dimensional  space.   Moreover,  our  method  gives  an  indication 
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about  stability  and  makes  it  possible  to  compare  different 
configurations  and  decide  which  one  is  the  most  stable. 

As  a  mathematical  model  of  the  Tokamak,  we  investigate 
the  magnetohydrodynamic  equilibrium  of  a  plasma  contained 
in  a  toroidal  region  D,  of  space  that  is  separated  by  a  sharp 
free  boundary  T      from  an  outer  vacuum  region   D~   bounded 
by  a  conductor  shell  C  .   The  plasma  is  assumed  to  be  a 
perfectly  conducting  compressible  fluid  at  constant  pres- 
sure  p   and  internal  energy   e  =  p/Cy-l). 

We  express  the  magnetic  field   B_   in  the  vacuum  region 
as  a  linear  combination 

B2  =  c,Vi[i,  +  C2^^2 

of  the  gradients  of  two  harmonic  functions  which  yield  unit 
currents,  and  the  magnetic  field  B,   in  the  plasma  as  a 
multiple 

B3  =  0374.3 

of  the  gradient  of  a  single  harmonic  function  with  unit 
current.   The  scalars   c,  ,02?  C3   are  determined  in 
such  a  way  that   B,  and  B-   yield  prescribed  fluxes. 

The  basic  idea  is  that  the  total  potential  energy 


E  =  I  j   B^  dV  +  ^  J   b2  dV  -f 

D2  D^  b^ 


e  dV 


is  stationary  for  an  equilibrium  position.  In  the  case  of 
a  stable  equilibrium,  the  energy  is  a  minimum  with  respect 
to  perturbations  of  the  free  boundary  T ,   but  a  maximum 
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with  respect  to  variations  of  ^^     ,    \p      ,    \p^    .      Thus  the 
problem  can  be  characterized  as  a  minimax  problem. 

To  determine  the  shape  of  the  free  boundary  T      we 
apply  the  idea  of  paths  of  steepest  descent  to  the  varia- 
tional principle,  using  the  requirement  that  the  first 
variation  of  the  energy 


=  j  (i  B^  -  i  bJ  -  p)  5„ 


(SE  =    (^  B^  -  =-  B,  -  p)  6n  dS 

r 

vanish  at  the  equilibrium  position.   We  thus  introduce 
an  artificial  time  parameter   t  ,   and  the  solution  is 
obtained  by  solving  a  time  dependent  system  of  equations 
which  converge  to  a  steady  state.   This  enables  us  to 
study  questions  of  equilibrium  and  stability  with  far  less 
computational  effort  than  would  be  necessary  if  we  examined 
dependence  on     physical  time  instead. 

We  have  developed  a  finite  difference  scheme  based 
entirely  on  a  discrete  version  of  the  variational  principle. 
We  establish  a  parametric  formulation  by  introducing  a 
fixed  cube  in  an  auxiliary  domain  [5] .   One  section  of  the 
cube  is  mapped  onto  the  vacuum,  and  the  remainder  is  mapped 
onto  the  plasma,  in  such  a  way  that  an  intermediate  plane 
is  mapped  onto  the  free  boundary  Y    .   We  seek  the  potentials 
^■1     /    ^2    '    ^"i      ^^  solutions  to  Neumann  problems  in  the  cube. 
The  difficulty  of  locating  the  free  boundary  is  overcome 
because  we  solve  all  equations  in  the  fixed  cube. 

The  coordinate  system  becomes  singular  along  a  closed 
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curve  inside  the  plasma  analogous  to  the  axis  in  cylindrical 
coordinates.   Difficulties  with  approximating  Laplace's 
equation  at  the  axial  curve,  with  Neumann  boundary  conditions 
and  with  avoiding  incompatibility  of  the  Neumann  problem  for 
^    ,^^,^^      are  all  solved  by  the  discrete  variational 
principle.   Furthermore,  the  method  converges  exclusively 
to  equilibria  that  are  physically  stable,  so  that  the 
convergence  of  the  process  is  a  measure  of  the  stability  of 
the  solution. 

The  fundamental  conclusion  is  that  the  method  can  be 
applied  successfully  to  fully  three-dimensional  configura- 
tions, and  that  it  yields  results  concerning  their  stability. 

The  computer  program  has  been  used  to  find  three- 
dimensional  magnetohydrodynamic  equilibria  both  with  and 
without  axial  symmetry.   We  have  found  a  great  number  of 
stable  configurations,  and  we  have  discovered  geometries 
without  axial  symmetry  that  are  shown  by  the  method  to  be 
more  stable  than  neighboring  axially  symmetric  cases. 

Even  though  we  have  considered  only  a  very  simple 
physical  model,  the  difficulties  arising  from  the  three  space 
space  dimensions  are  enough  to  nearly  exhaust  the  capacity 
of  the  CDC  6  600  computer. 
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2 .   Magnetohydrodynamics  Model. 

We  are  interested  in  the  problem  of  containment  of 
an  ionized  gas,  or  plasma,  by  a  strong  magnetic  field. 
Therefore  our  aim  is  to  investigate  the  equilibrium  of  a 
conducting  fluid  in  the  presence  of  a  magnetic  field. 

Magnetohydrodynamics  is  one  of  the  models  describing 
this  behavior.   The  model  is  appropriate  because  it 
combines  reasonable  physical  accuracy  with  mathematical 
simplicity. 

In  this  approach,  the  plasma  is  considered  to  be  a 
classical  perfectly  conducting  fluid,  characterized  by  its 
density   p  ,   its  entropy  S   and  its  flow  velocity   u. 
The  pressure  is  then  a  function  of   p  and  S.   The  magnetic 
field   B   and  the  current   J   interact  with  the  plasma  and 
exert  a  Lorentz  force   J  x  b   on  the  plasma.   In  order  to 
obtain  a  simple  mathematical  model,  we  neglect  displacement 
currents,  electrostatic  forces  and  charge  density.  Similarly, 
the  effects  of  viscosity  and  heat  conduction  will  be  omitted 
[1] .   The  basic  principles  which  describe  magnetohydrodynamic 
flow  are  Maxwell's  equations  and  Newton's  second  law  of 
motion,  together  with  the  laws  of  conservation  of  mass 
and  energy. 


-5- 


2.1.   The  Lundquist  Equations 

According  to  Ohm's  Law  for  a  perfect  conductor  we  have 

(1)  E  =  B  X  u 

where   E   is  the  electric  field  vector.   Since  the  displace- 
ment current  is  neglected,  from  Maxwell's  equations  we  can 
express   J   in  terms  of  the  magnetic  field  by 

(2)  J  =  V  X  B 

and  the  Lorentz  force  per  unit  volume  of  fluid  becomes 

J  X  B  ^  -  Bx  (VxB)  . 

It  is  now  possible  to  eliminate  all  of  the  electro- 
magnetic variables  except  B  .  From  Maxwell's  equations 
we  have 

(3)  V-B  =  0 

and  Faraday's  law  of  induction  gives 

1^  +  VxE  =  0 

Using  (1),  (2)  and  putting   d/dx  =  3/9t  +  u*V,   Faraday's 
law  reduces  to 

(4)  ^  +  (V-u)B  -  (B-V)u  =  0 
Newton's  second  law  of  motion  yields 

(5)  p  ^  +  Vp  +  Bx(VxB)  =  0 

dx 

where  the  last  term  represents  the  Lorentz  force  exerted  on 
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the  fluid  by  the  magnetic  field.   The  continuity  equation 
for  the  fluid  is 

(6)  k^^   ^.^  =    0 

p  dT 

We  assume  that  the  pressure  and  density  of  the  fluid 
are  connected  with  the  entropy  by  an  equation  of  state  of 
the  form 

(7)  p  =  A(S)p^ 

where  y      is  the  ratio  of  specific  heats.   The  law  of 
conservation  of  energy  then  states 

(8,  If  =  0 

The  system  of  equations  (4) -(8)  together  with  the 
initial  condition  (3)   are  known  as  the  Lundquist  equations. 

2 , 2    Equilibrium  Equations 

In  this  paper  we  shall  be  concerned  with  determining 
configurations  of  static  equilibrium  in  magnetohydrodynamics 
for  a  perfectly  conducting  fluid  which  is  contained  by  a 
magnetic  field.   In  this  case  there  is  no  flow  and,  there- 
fore ,    u  =  0  . 

As  a  result  of  Ohm's  law,  the  electric  field  E  must 
also  be  zero.   Next  equation  (6)  states  that 

and  p   must  be  a  constant.   Thus,  for  equilibrium,  the 
Lundquist  equations  reduce  to  the  system 
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(9)  V-B  =  0 

(10)  Vp  +  Bx (VxB)  =  0 

In  order  to  derive  boundary  conditions  to  be  imposed 

at  the  interface  between  two  distinct  media,  we  first  have 

to  convert  (9)  and  (10)  into  integral  relations  which  will 

determine  the  jumps  of  B  and  p   across  such  a  surface  of 
discontinuity. 

Let   D   be  an  arbitrary  region.   Applying  the  divergence 
theorem  to  (9)  we  have 

(11)  B-n  dS  =  0 

9D 

where   9D   is  the  boundary  of  D,   n   its  normal  and   dS   the 

area  element  on   3D. 

A  combination  of  (9)  and  (10),  and  another  application 
of  the  divergence  theorem  yields 


(12) 


; 


[ (p  +  y  B^)n  -  (B-n)B]  dS  =  0 


3D 
Assume   D   lies  across  an  interface   F   which  separates 

two  distinct  media.  Both  of  the  integral  relations  (11) 
and  (12)  remain  valid,  even  if  the  domain  taken  is  arbi- 
trarily small.   Under  the  assuitiption  that   B   is  bounded, 

2 
we  obtain  that   B*n  =  0  on  r  and   B  /2  +  p   must  be 

continuous  across   r  . 
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2 . 3    Description  of  the  Physical  Problem 

We  shall  describe  the  equilibrium  problem  which  will  be 
considered.   Let   D,   be  the  region  occupied  by  the  plasma, 
surrounded  by  a  vacuum  region   D^  .   The  common  boundary 
separating  the  plasma  and  vacuum  regions  will  be  denoted 
by   r  ,   and  the  outer  boundary  of   D„   will  be  denoted  by  C. 
We  shall  assume  further  that  the  location  of   C   is  fixed, 
and  that  it  has  a  toroidal  shape.   The  location  of   F   is 
determined  by  the  equilibrium  position  of  the  plasma 
(Figure  1) . 

We  shall  make  the  extra  assumption  that  inside  the 
plasma  the  current  is  zero;  i.e.,  all  the  current  is 
carried  on  the  interface  T      (skin  current) .  This  assump- 
tion will  be  justified  later  in  Chapter  3,  by  an  applica- 
tion of  the  variational  principle. 

It  follows  from  (2),  (9)  and  (10)  that  B^^  ,  the 
magnetic  field  inside  the  plasma  satisfies 

(13)  V-B^  =  0 

(14)  VxB^   =    0 

and  that  the  pressure  p  reduces  to  a  constant. 

No  current  exists  in  the  vacuum.  Hence,  it  follows 

from  (2)  and  (9)  that  B-  ,   the  magnetic  field  in  the 
vacuum,  satisfies" 

(15)  V-B2  =  0 

(16)  VXB2  =  0  . 
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At  the  common  boundary  T      between   D,  and  D„  ,   the 
jump  conditions  which  we  derived  in  the  last  section  yield 
the  boundary  conditions  on   F 

12        12 

(17)  i  B^  +  p  =  i  B2 

(18)  Bj^-n  =  B2*n  =  0 

where   n   is  the  normal  to  T    ,   and   •   denotes  the  scalar 
product.    Similarly,  at  the  outer  boundary  C   of  the 
domain   D„   we  have  the  boundary  condition 

(19)  B2-n  =  0 

where   n   is  the  normal  to  C  . 

In  addition  to  the  differential  equations  (13) -(16), 
and  the  boundary  conditions  (17) -(19),  the  free  boundary 
problem  requires  the  specification  of  an  appropriate  number 
of  periods.  The  exact  number  is  determined  by  the  geometry 
of  the  problem,  and  in  the  case  of  a  torus  is  three.   These 
periods  may  be  specified  in  two  ways,  either  as  fluxes  or 
currents. 

Let  fi,  ,    0,2    ,    f^o   be  the  surfaces  shown  in  Figure  1. 
Then  the  corresponding  fluxes  are  defined  by 

f 


(20)  F.  = 


B'dS  ,      i  =  1,2,3 


1 


where  dS    is  the  vector  area  element  on  9..     .      Let 
C,,C„,C-   be  the  closed  loops  shown  in  Figure  1.   Then 
the  corresponding  currents  are  defined  by 
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(21)  I.  =  c)   B-d£  ,      i  =  1,2,3 

C. 

1 

where  di      is  the  vector  length  element  on  C .  . 

Since   B,  and  B^   are  curl  free,  we  can  find  scalar 
potentials   a   and   3   such  that 

B,  =  Va 

(22)  ^ 

B2  =  VB 

It  follows  that  equations  (13) -(16)  reduce  to 

Aa  =  0   in   Di  / 

(23)  ^ 

A3  =  0   in   D, 


2 

and  the  boundary  conditions  (17) -(19)  reduce  to 


(24) 


dn 
dn 


on  r, 

(25)  1^  =  0 


8n 


on  C ,  and 


(26)  I  (Va)^  +  P  =  i  (^2)^ 

on  r.   Given  the  position  of  the  boundary,  the  differential 
equations  (23)  and  the  boundary  conditions  (24) ,  (25)  toge- 
ther with  the  specification  of  the  three  periods  (either  (20) 
or  (21))  determine  the  solutions   a   and   6   uniquely  up 
to  a  constant. 
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The  equilibrium  condition  (26)   represents  a  balance 
of  forces  on  the  boundary  T    ,      and  it  determines  the 
equilibrium  position  of  the  free  boundary  r  . 
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3.   Variational  Principle  for  the  Boundary 

Of  even  greater  importance  than  the  computation  of 
configurations  of  magnetohydrodynamic  equilibrium  is  the 
determination  of  their  stability.   We  choose  a  method  which 
gives  some  indication  of  stability,  for  the  numerical  scheme 
we  use  is  based  on  the  idea  of  making  the  energy  a  minimum 
by  the  method  of  steepest  descent. 

As  described  in  the  previous  chapter,  to  determine  the 
solution  of  the  differential  equations  we  must  prescribe 
three  periods,  which  may  be  specified  as  fluxes  or  currents. 
However,  Maxwell's  equations  with  boundary  conditions 

B'n  =  0  ,    (E+uxB) •£  =  0 

at  the  interface  moving  with  fluid  velocity  u  (where  n  and 
£   are  unit  vectors  in  the  normal  and  tangential  directions 
respectively)  imply  that  the  fluxes  are  constants  of  the 
motion  [8].   To  prove  this  statement  let  Q. .   (t)   be  the 
position  of   $1.   at  time   t,   and  let  F.  (t)   be  the  flux 

across  Q.      at  time  t.  Then 

1 

dF. 


dt 


d 

dt        J 

B'dS 

^^(t) 

n,U) 

If  +     (V.B)u 

•dS  + 


Bxu'dJl 


L^(t) 


where   u  is  the  velocity  of  the  surface  Q.    ,  dS  is  the 
vector  area  element  on  Q.     ,      L.   is  the  closed  loop  which 
bounds  Q.     ,  and  d£  is  the  directed  length  element  on  L.. 
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From  Maxwell's  equations,  and  from  the  initial  condition 
V'B  =  0,  we  have  that   V'B  =  0   for  all  t.   Using  Faraday's 
law  of  induction,  the  last  equation  reduces  to 


dp. 

^^  =   -         I         VxE'dS   +       I  Bxu'dil 

f^j_(t)  L^(t) 


(E+uxB) 'dl      =      0 


L^(t) 


Thus  in  order  to  obtain  a  physically  meaningful  equili- 
brium, we  must  choose  a  formulation  in  which  the  fluxes  are 
kept  fixed  for  different  positions  of  the  free  boundary  r. 

If  we  assume  that  the  problem  has  a  certain  type  of 
symmetry,  for  example,  axial  or  helical  symmetry,  we  can 
introduce  a  flux  function  which  satisfies  Laplace's  equation 
with  given  boundary  values  (the  prescribed  fluxes) .   In  the 
general  case   of  a  fully  three  dimensional  problem,  no  such 
single  flux  function  exists  and  we  are  forced  to  work  with 
the  scalar  potentials  a  and  B  . 

In  this  chapter,  a  variational  principle  is  formulated 
to  find  the  equilibrium  position  of  the  plasma  boundary  which 
corresponds  to  a  minimum  of  energy  of  the  system. 

3.1.   The  First  Variation  of  the  Energy. 

We  will  determine  the  position  of  the  free  boundary  by 
applying  the  method  of  steepest  descent  to  the  following 
variational  principle   associated  with  the  free  boundary 
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* 

problem.   Among  all  admissible  boundaries  T      ,    each  of 

which  defines  domains  D,  and  D   in  which  we  can  find 
potentials  a  and  6  satisfying  the  equations  (23)  and  the 
boundary  conditions  (24) ,  (25)  subject  to  the  constraints 
(20) ,  the  appropriate  free  boundary  T      is  characterized 
by  the  property  that  the  energy  E   of  the  system  is  a 
minimum  for  a  fixed  mass   M  of  plasma. 

The  appropriate  expression  for  the  energy  of  the 
system  is  [2] 


(27)       E  =  J 


2       1 
B^  dV  +  ±- 


B^  dV  +  I  e  dV 


^1  °2         °1 


where   e  =  p/(Y-l)   is  the  internal  energy  per  unit  volume 
of  the  fluid.   The  first  term  represents   the  magnetic 
energy  E^  in  the  plasma,  the  second  term  represents  the 
magnetic  energy  E_  in  the  vacuum,  and  the  last  term  E^ 
represents  the  internal  energy  of  the  fluid. 

To  each  position  of  the  free  boundary  r  we  have 
associated  a  unique  value  of  the  energy   E(r)   by  requiring 
that  the  magnetic  fields  satisfy 

VxB.  =  0 

^  i  =  1,2 

V'B.  =  0 

1 

subject  to  the  constraints  (20)  ,  boundary  conditions  (18)  , 
and  that  the  pressure  p  =  constant.   For  this  procedure  to 
be  physically  meaningful,   E(r)   must  correspond  to  the 
minimum  energy  for  all  admissible  states.   To  prove  this 
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let   E(B.,p)   be  the  energy  corresponding  to  a  state 
characterized  by  magnetic  fields   B.  and  pressure  p, 
subject  to   V»B.  =  0  ,      p  dV  =  M,   constraints  (20)  and 
boundary  conditions  (18). 

If  we  minimize   E(B.,p)   within  this  class  of  states, 
the  Euler  equations  for  the  variational  problem  yield 

VxB^   =0         i  =  1,2 
p  =  constant 

so  that 

E  (r)  =  min  E(B.  ,p)  . 
B^,p    ^ 

For  computational  purposes,  this  variational  principle 

is  not  practical  in  the  case  of  a  fully  three  dimensional 

problem.   Instead,  we  use  the  reciprocal  variational 

principle  in  which  we  maximize  the  magnetic  energy  subject 

to  the  constraint   VxB.  =  0.   That  is,  we  use  the  same 

1 

boundary  conditions,  constant  fluxes  but  we  replace  the 
constraint   V*B.  =  0  by  the  constraint   VxB.  =0.   In  this 
case,  the  Euler  equations  yield 

V«B.  =  0  ,        i  =  1,2 
1 

so  that  the  maximum  is  attained  by  an  admissible  state,  and 

because   VxB.  -    0   the  corresponding  energy  is  minimum  among 

all  admissible  states.   In  other  words,  both  procedures  yield 

the  same  magnetic  fields  corresponding  to  the  minimum  of 

energy  E(r)  among  all  admissible  states. 
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This  is  why  we  are  justified  in  assuming  that  the 
current  inside  the  plasma  is  zero,  and  we  can  express  B 
as  the  gradient  of  a  scalar  function. 

Going  back  to  the  variational  principle  for  the  energy 
E(r)   as  a  functional  of  the  boundary  position  we  have 

Theorem.   Let   5v   be  an  arbitrary  perturbation  of 
the  free  boundary  T      along  its  outer  normal  (cf .  Figure  2) . 
Then 

'  '1  „2    1  „2 


(28)         6E  = 


^J  ^2    ~   J  K    ~   P^  ^^    '^^ 


r 

where   dS   is  the  area  element  on  T    . 

Proof:   Let  ^,      be  the  potential  for  the  vacuum  region 
corresponding  to  the  currents   I,  =1,   l2=0,   i.e., 
4^,   satisfies  Laplace's  equatioi 
conditions,  and  its  periods  are 


4^,   satisfies  Laplace's  equation,     Neumann  boundary 


I,  =  1  Vi|y^'d£  =  1  ;     ±2    =         Vi|;^«d£  =  0  . 
^1  ^2 

Similarly,  let   4^2  ^^   ^^®  potential  for  the  vacuum  region 
corresponding  to  I,  =  0  ,   I^  =  1. 

It  is  possible  to  express  the  potential   3   for  the 
vacuum  region  in  terms  of   ij^,  and  i/i-  ^.s    follows: 

(29)  B  =  c^ip^   +  C2ip2 

where   c,  and  c-   are  scale  factors  (physically  they 
correspond  to  the  currents  for  the  potential  3)   which 
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will   be   determined     so   that      B      satisfies    (20)    for   i   =   1,2 
We    define    the    inductance   matrix      A     by 


(30)  A  -    a^j    =      J    ViPj-dS    ,  i   =    1,2 

n. 

J. 

where   dS   is  the  vector  area  element  on   f^ .  .   From  (29) 
and  (30)  we  have 

(31)  f  =  Ac 

where   f  is  a  vector  whose  components  are  the  prescribed 
fluxes   F, ,F_   and   c  is  a  vector  whose  components  are  the 
scale  factors   c, ,c_. 

An  application  of  Green's  theorem  shows  that 


(32)  a^j  =   J  V.j;^-Vii;^  dV  ,       i  =  1,2;  j  =  1,2 

^2 
and  furthermore  from  (2  9)  and  (32)  we  have 

(33)  I  (V6)^dV  =  c'Ac  >  0 

^2 
where   c'  is  the  transpose  of  c. 

These  last  two  equations  prove  the  well  known  result 
that  the  inductance  matrix  A  is  symmetric  and  positive 
definite. 

We  can  compute  the  first  variation  of  the  magnetic 

energy 

E2  =  I  J  (VB)^dV  =  I  c'Ac  =  I  c'f  =  i  f'c  . 

°2 
Since   f  is  constant,  (31)   yields 
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and 


so   that 


6E2   =   J   f  (Sc   =   I  c'A'6c   =   J  c'A6c 


6c  =   -   A      6a  c 


(34)  6E2    =    -   ^  c'6A   c      . 


The   expression    for      5 A  is   given  by   the   following 
Leinma.      We  have 


6a^  .    =    -         (Vi(j^«Vi|j  .)    6v    dS 

r 

where   dS   is  the  surface  area  element  on  r  . 

Proof  of  the  lemma:   We  will  define  a  variation  of  D- 

* 
onto   Dp   by  introducing  the  notion  of  an  interior  variation 


of   D2   [Ij.   Let 


S .  =  S  (x,y,z)  ,     j  =  1,2,3 


be  functions  possessing  continuous  partial  derivatives  of 
a  sufficiently  high  order  in  some  neighborhood  of  D2  • 

For  small  enough  choices  of   e  ,  the  transformation 

* 
X   =  X  +  e  S,  (x,y ,  z) 

(35)  y   =  y  +  e  S2  (x,y,  z) 

* 
z   =  z  +  e  S,  (x,y ,  z) 

maps   D  onto  D  in  a  one-to-one  fashion.  The  location  of 

the  fixed  boundary  C  of  the  domain   D_   must  be  identical 

*     * 
to   that  of  the  fixed  boundary  C  of  D_  .   Hence, 

Sj^(x,y,z)  =  S2(x,y,z)  =  S2(x,y,z)  =  0 
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for  all  points   (x,y,z)   lying  on  the  boundary  C  .   Let   6v 
be  the  displacement  of  the  free  boundary   T   under  such  a 
transformation.   Then  from  (35)  we  have 

(36)  6v  =  e(S^  If  +  S2  U+  S3  If)  . 

*  "  * 

We  introduce  on   D2   potentials  \l).      which  satisfy 

Laplace's  equation  with  Neumann  boundary  conditions  and 
with  periods  identical  to  those  of  it.    ,    i.e. 

*  * 

Alp.    =    0      in  ^2    '  1   =    1,2 

8ijj .  *         * 

^-i  =OonC,r,  i=l,2 

dn 


J    a. 


•dJl   =    I     Vij;^-d£    ,  i    =    1,2;     j  =  l,2 


C*  C. 

D  J 


and  we   define 


* 


ViJ; .  •  VtiJ  .    dV      ,  i    =    1,2;    j  =  l,2 


* 


°2 


the    inductance   matrix   in      D2    ,    so    that 


6a. .  =  e  lim 


* 
a.  .-a.  . 


^^       £-0 


** 


We  introduce  comparison  functions  ^^      defined  on  D2  by 

ie  ic         ic         ic         "k 

^^    (x  ,y  ,z  )  =  iiJ^(x,y,z)  ,      i  =  1,2 


and  let 


a.  .  =    Vi|;.  'Vilj.   ,  1  =  1,2  ]  =  1,2 


°2 
Direct  calculation  shows  that 
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C.  C. 


Vi/j.  •d£    +    O(e^) 


**  **  * 

so    that        \l>.      and      ii.     ,    and  consequently      ^.      and   ip .     ,    have 

the  same  periods. 

**  * 

The  function  ^.         is  not  harmonic  in   D_  ,  but  we 

* 
can  appeal  to  Dirichlet's  principle,  which  states  that  a.  . 

is  stationary  for  solutions  of  Laplace's  equation  with 

Neumann  boundary  conditions.   It  follows  that 

**     *        2 
a.  .  =  a.  .  +  0(e) 

so  suffices  to  prove  that 

a**  =  a.  .  -    (Vi|;.  -Vii).)  6v  dS  . 

r 

This  last  result  follows  from  a  direct  but  lengthy  calcula- 
tion which  uses  the  fact  that  i> .  ,'^  .      are  harmonic  functions 
and  satisfy  Neumann  boundary  conditions  [1] .   This  completes 
the  proof  of  the  lemma. 

Continuing  the  proof  of  the  main  theorem,  we  have  from 
(34)  and  the  lemma 

6E2  =  J  I  (c^ViJ;^+2c^C2ViJj^«Vi|;2+C2Vt|j2)  5v  dS 

r 

which  by  (29)  reduces  to 
(37)      "^^2  "  i  J  ^^^^^  '^^  ^^  =  i  J  ^2  "^^  ^^  • 

r  r 

We  follow  the  same  procedure  to  compute  the  first 
variation  of  E,  ,  the  magnetic  energy  in  the  plasma.   In 
this  case  only  one  period  is  needed  to  determine  the 
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potential   a  ,   and  the  inductance  matrix  reduces  to  a 
scalar.   Let  ^^      be  the  potential  for  the  plasma  region 
corresponding  to   I^  =  1.   Then,  instead  of  (29) -(32)  we 


have 


(38)  a  =    C3IP3 


(39)  333   =        ViJj^-dS 

(40)  F3  =  33303       ^:'  " 

(41)  333  =  I  (ViP3)^  dV  . 

We  compute  the  first  variation  of 

^1  =  H  (^«)'  d^  =  J  ^3^3 
°i 

From  (39) ,  and  the  fact  that   F-  is  constant,  we  have 


'^^l  =  "  2  ^3  "^^33 


Applying  the  lemma  to  D,  we  have 


6333  =  I  (ViJ;3)^  6v  dS 

r 

where  the  positive  sign  is  due  to  the  fact  that  6v  is 
positive  along  the  outer  normal  to  the  domain  D,  ,  but 
6v  is  positive  along  the  inner  normal  to  the  domain  D^ 


And  finally 


(42)    6e,  =  - 


1  T  ,r,  ^2  ,    ,„       1  f  „2 


(Va)   5v  dS  =  -  ^   B"  6v  dS 


1      2  J  ^'"'  ^.    ^^   -        2    j    "1 
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Next  we  compute  the  first  variation  of  the  internal 
energy  E^  ,   for  a  fixed  mass   M  of  plasma,  with  equation 
of  state  (7) .   We  have 

°1 


(43)  p  =   k  p 


Y 


p  dV   =   M 
°1 
where  p  and  p   are  constant  throughout  the  plasma.  From 
(43)  we  find 

6E3  =  ;^  J  p  6v  dS  +  ;j33-  j   6p  dV 


r      .        D, 


6p  =  I^^ 
^      P 


6M  =    6p  dV  +    p  6v  dS  =  0  . 

D,      r 

since   p  and  p  are  constant  in  D,  ,  the  last  equation 

implies 

I  p  6v  dS  +  I   2_A^  dV  =  0 

'  °1 

which  can  be  written  as 

.  6p  dV  =  -  Y    P  <5v  dS 
so  that         1 

6E3  =  ;p3-  J  p  6V  dS  -  ;j^  I  p  6V  dS 
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(44)  SE^    =  -    P  6v  dS 

r 

finally  combining  (37) ,  (42)  and  (44)  we  have 

6E  =  6E^+6E2+6E3  =  j  (|  B^  -  j  bJ  -  p) 6v  dS  . 

r 

3.2.   The  Method  of  Steepest  Descent. 

To  have  equilibrium,  it  is  a  necessary  and  sufficient 
condition  that  the  first  variation  &E      must  be  equal  to  zero. 
Since  6v  is  arbitrary,  the  integrand  in  (28)  must  be 
identically  zero.   This  is,  of  course,  equivalent  to  the 
free  boundary  condition  given  by  (26)  .  The  expression  for 
the  first  .variation  provides  us  with  a  method  of  updating 
the  position  of  the  free  boundary  T    so  as  to  satisfy  (26) 
and  find  a  minimum  for  E. 

We  choose  a  path  of  steepest  descent  by  setting 

(45)  6v  =  -  X(|  B^  -  I  bJ  -  p) 

where   X  >  0   represents  a  convergence  factor.  Then   6E  will 
be  negative,  each  new  position  of  the  boundary  T      will 
correspond  to  a  lower  energy,  and  by  interation  we  will  reach 
a  minimum  for  which   6E  =  0.   If  such  a  minimum  does  not 
exist,   for  example  if  the  equilibri\im  is  unstable,  the 
process  will  not  converge. 
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3.3.   Characterization  of  the  Equilibrium  Position  as  a 
Saddle  Point  for  the  Energy. 

We  have  discussed  the   behavior  of  the  energy  as  we 
vary  the  position  of  the  free  boundary.   We  shall  now 
describe  the  behavior   of  the  energy  when  we  hold  the 
free  boundary  fixed  but  vary  the  values  of  the  potentials, 

We  define  on  D„ 

^       * 


A  =  {^1  I  I^  =  0,  I^  =  1} 
>  =  (*2  1  ^1  =  1'  ^2  =  0> 


and  on  D, 


Let 


±[   =    ^^l    \    h=    ^^ 


A*  =    Vi/j*-V4j*  dV,        i=l,2;j  =  l,2 


1    D 


and  °2 


^33  =  J  ^^3-^'^3 

°1 

* 

Then  the  corresponding  expression  for  the  energy  E   as  a 

* 
functional  of  i) .      is 

E  {i>^)    =   i  [f  (A  )  -^f  +  F3(a33)  ■^F3]  +  ;pi  V  , 

where  V  is  the  volume  of  D,  .   As  we  vary   lJ;'^   in    J~  •     , 

*  * 

E   is  maximum  when  i) .      is  a  solution  of  Laplace's  equation 

with  Neumann  boundary  conditions.   We  can  illustrate  this 
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result  by  examining  the  axially  symmetric  case,  for  which 
the  matrix  A  is  diagonal.   In  that  case,  the  expression  for 
the  energy  reduces  to 


^3 


*     1   ?•  -1       P 

(ip.)  =  4   I   F.  (a.  .)  -^F.  +  -~ 
1      2  .^^   1   11      1    Y-l 


.  V 
i=l  -      -  -    Y- 


imum 


and  by  Dirichlet's  principle  we  have  that   a. .   is  mini 

for   ij; .   solving  Laplace's  equation  with  Neumann  boundary 

* 
conditions.   Therefore,  the  last  equation  shows  that  E 

is  maximum. 

Let   E(r)   denote  the   value  of  this  maximum  for  a  given 

position  of  the  boundary  T    .      As  we  vary  T,      a  stationary 

point  of   E(r)   corresponds  to  an  equilibrium  position  and 

a  minimum  of  E(r)  corrresponds  to  a  stable  equilibrium.  If 

such  a  minimum  exists,  the  energy  at  the  equilibrium 

position  is  given  by 

*   * 

E  =  min  E(r)  =  min  max  E  (ip.) 

r        T      ip*  ^ 

1 
It  follows  that  the  problem  viewed  as  a  whole   is  a 

minimax  problem.   In  other  words,  the  equilibrium  position 
corresponds  to   a  saddle  point  of  the  energy.   This  charac- 
teristic of  the  solution  is  due  to  the  fact  that  we  have 
used  a  formulation  in  terms  of  a  potential   instead  of 
a  flux  function. 

3.4.   Remarks  on  the  Method. 

The  advantages  of  this  method  are  the  following:  We 
have   reduced  the  problem  to  the  solution  of  three  scalar 
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equations.   This  formulation  allows  us  to  express  the 
first  variation  of  the  energy  in  terms  of  the  first  variation 
of  Dirichlet's  integral,  so  that  we  can  use  Dirichlet's 
principle.   Furthermore,  from  the  computational  point  of 
view,  we  have  to  solve  Laplace's  equation   with  fixed 
current  ,  which  are  the  natural  periods  for  the  potentials. 

The  disadvantage  is  that  we  have  a  minimax  problem, 
and  this  requires  that  for  each  position  of  the  boundary 
Laplace's  equation  must  be  solved  very  accurately. 

The  solution  is  obtained  by  the  following  steps: 

1.  We  guess  a  position  of  the  free  boundary,  and  solve 
for  the  potentials  ^i,  ,il)~,\p^. 

2.  We  compute  the  inductance  matrix  and  use  it  to  compute 
the  scale  factors   c,,c  ,c   so  that  (31)  and  (40) 

are  satisfied. 

3.  We  use  the  method  of  steepest  descent  to  find  a  new 
position  of  the  boundary  T    .   The  process  is 
terminated  when  the  free  boundary  condition  (17)  is 
satisfied  within  a  given  accuracy. 

The  next  chapter  deals  with  a  numerical  method  for 
computing  the  solutions  for  the  potentials  4^,  ,  ijj„ ,  ijj-,  for  any 
given  boundary   r   and  a  numerical  method  for  updating 
the  boundary  position. 
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4 .   Discrete  Variational  Principle  for  Laplace's  Equation 

The  basic  difficulty  in  developing  a  numerical  scheme 
for  a  free  boundary  problem  is  to  be  able  to  handle  the 
changing  shape  and  location  of  the  free  boundary.  If  the 
problem  involves  two  independent  variables,  conformal 
mapping  techniques  can  be  used  [3,4]  to  solve  the  problem 
in  a  fixed  auxiliary  domain  R,  and  then  the  region  of 
solution  is  given  as  the  conformal  image  of  R.  Thus,  the 
problem  of  computing  the  location  of  the  free  boundary  is 
transferred  to  that  of  finding  such  a  conformal  mapping. 

In  the  general  case  of  three  independent  variables,  no 
such  conformal  mapping  exists,  but  the  basic  idea  of  mapping 
the  flow   region  into  a  fixed  auxiliary  domain  R  can  be 
used.   In  this  case,  since  the  mapping  is  not  conformal, 
the  equation  in  R  will  be  more  complicated,  which  in  turn 
means  that  finding  the  solution  will  require  a  greater 
amount   of  computation.   Nevertheless,  our  experience  shows 
that  solving  the  problem  directly  in  the  flow  region  is 
practically  impossible  due  to  the  changing  location  of  the 
free  boundary. 

The  first  step  is  to  choose  an  appropriate  coordinate 
system,  so  that  the  mapping  and  the  resulting  numerical 
scheme  in  the  auxiliary  domain  are  as  simple  as  possible. 
Thus  we  keep  the  amount  of  computation  to  a  minimum. 

The  second  step  is  to  write  a  numerical  scheme  for  the 
partial  differential  equation   and  boundary  conditions  in 
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the  auxiliary  domain.   There  are  three  problems  to  be 
considered  here.   The  first  one  is  how  to  write  the  equation 
on  the  boundary   using  Neumann  boundary  conditions.   The 
second   is  how  to  get  a  compatible  system  of  difference 
equations.   If  we  write  the  difference  equations  as  a 
system  of  the  form  Ax  =  b,  the  coefficient  matrix  A  has 
zero  determinant  because  a  constant  is  a  solution  of  the 
homogeneous  problem.   Therefore,  in  order  for  a  solution 
X  to  exist,  the  system  Ax  =  b  must  be  compatible.  The 
third  problem  is  how  to  write  the  difference  equations  for 
singular  points  of  the  coordinate  system  like  those  in  the 
case  of  cylindrical  coordinates  where  r  =  0  . 

These  problems  can  be  solved  by  using  a  discrete 
■  ariational  principle  for  Laplace's  equation.   Since 
Neumann  boundary  conditions  are  the  natural  boundary 
conditions  for  Laplace's  equation,  we  can  get  the  equations 
for  boundary  points  by  varying  the  values  of  the  solution 
on  the  boundary,  and  they  will  automatically  satisfy  the 
boundary  conditions.   Since  the  difference  equations  are 
derived  from  a  variational  principle,  the  resulting  system 
will  be  compatible.   The  variation  of  the  solution  at  the 
singular  points  yields  the  appropriate  mean  value  theorem 
for  the  solution  at  the  singular  points. 

The  third  step  is  to  write  a  numerical  scheme  for  the 
updating  of  the  position  of  the  free  boundary.   Again,  the 
appropriate  approximation  is  suggested  by  the  variational 
principle. 
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4.1.   Coordinate  System. 

We  choose  an  orthogonal  system  of   coordinates  r,()),9 
as  follows  (cf.  Figure  3):  First,  we  draw  a  circle  of  radius 
a   corresponding  to  the  major  radius  of  a  torus,  which  we 
center  at  the  origin  on  the  plane   z  =  0.   Let  the  angle 

8  be  the  same  as  in  cylindrical  coordinates.   In  the  plane 

9  =  constant,  with  origin  at  the  intersection  with  the  circle 
of  radius  a,  we  use  a  polar  coordinate  system  r,  4). 

We  assume  that  the  circle  of  radius   a   on  the  plane 
z  =  0   is  always  inside  the  plasma  region  D,  .   This  restric- 
tion can  later  be  lifted  if  we  replace  the  circle  of  radius 
a  by  any  given  closed  curve  in  space,  and  we  define  a  polar 
coordinate  system  on  the  plane    6  =  constant  with  origin 
at  the  intersection  with  the  given  curve.   In  this  case, 
the  restriction  would  be  that  the  given  curve  has  to  be 
inside  the  plasma  region. 

The  transformation  of  coordinates  we  have  described  is 
given  by 

X  =  (a  +  r  cos  (}))  cos  9 
y  =  (a  +  r  sin  (f)  sin  9 
z  =  r  sin  6 

and  its  Jacobian   J  is 

The  transformation  is  one-to-one  for  0  <  r  <  a.  We  have 
singular  points  for  r  =  0  (in  our  case  r  is  always  less  than  a) 
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The   expression    for    the   gradient   is    given   by 
Vij;(r,(}),8)    =    ^   e^   +  -^  e^   + 


r   r        r        <^        a+r   cos    (fi      6 

/N  /V  y\ 

where   ^r  '  ^4  '  ^6  ^^^   unit  vectors  in  the  r,4>,9  directions, 

A  general  surface  such  as  the  free  boundary  will  be 
given  by  r  =  f  ((t),9)  ,   and  f  =  constant   corresponds  to  an 
axially  symmetric  torus  with  circular  cross  section. 

4.2.   Variational  Principle  for  Laplace's  Equation  in  the 
Vacuum  Region. 

Let  the  equations  for  the  outer  boundary  C  and  the  free 
boundary  T      be  given  by: 

C:   r  =  x((}),  0) 

r:   r  =  f((t.,0) 

* 

We  define  a  mapping  from  the  vacuum  region  D^    onto  D_  by 


s  = 


r-f((}),9) 


x((t),e)  -  f((t),e) 

(46) 

9=6 
for        f£r£x;       Ol't'l  2it;       0   _<    9    <_  27r;    so   that 

D^    =    {(3,4), 9)     I     O^s^l,    0    <    (i)    <_  2iT,    0    <_   9    £   2-n} 

and  the  corresponding  images  of  C  and  V    are 

* 

C      =    {(s,(}),6),|    s=l;    0<_(i)<_2TT;    0<_9<_2Tr} 

r      =    {(s,<j),9)     I    s    =    0;    0    <    4)    <    27r;    0    <    9    <    2tt}    . 
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In  other  words,  D^  is  mapped  onto  a  rectangle  D„  ,  and 
the  surfaces  C  and  T  are  mapped  onto  planes  s  =  1  and  s  =  0 
(cf .    Figure    4) . 

Let    i|^  (r  (s,  4) ,  0)  ,  (J),  9)    =   g{s,(|),6);      we  want   to    solve 

Laplace's   equation    in    D„    ,    with  Neumann   boundary    conditions 

*    * 
on  C  ,  r   and  periodicity  conditions 

^\li'dl   =   g(s,(t)+2TT,  e)  -  g(s,(t),e)  =  I 

(47)       ^ 


ViJ;-d£   =    g(s,(}),  e+27r)    -    g{s,(p,Q)    =   !„ 
^2 
We   now   state    the    following   variational  principle    for  Laplace's 
equation.      Let     J~    be    the    class    of   all    functions    tp   with 
continuous    second   derivatives    in   D^    /    continuous    first 
derivatives   on   the   boundary,    and  periodicity   conditions 


^1  S 

Let  iin  ^    J"     be  the  function  that  renders  the  Dirichlet 
integral, 

F{i|^)  =  J  (V4;)^  dV 

°2 

stationary.      Then      i)^      satisfies   Laplace's   equation   in  D-    , 

and    9i|jQ/3n   =    0      on   the   boundary,    for 

6F(i|Jq)    =    -         Ai|Jq    6i|j    dV  +  j^   6(|j    dS   =    0 

D2  9D2 

and   since    6\J^    is    arbitrary,    we  have    Aijj-   =    0    in   D-    /    and 
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dipQ/dn   =    0    on    dD^  . 

We   now   use   this   variational  principle    to  derive   the 

* 

equation  for  g  in  D- •   Using  (46)  we  write  Dirichlet's 

integral  in  Dj    in  the  form 


F(g)  = 


(Vi|;)   dV 


(48) 


where 


D, 


=  J  (agg+bg^+cgQ+2dg^g^+2egggQ)  ds  d(|)  de 


r(s,cj),0) 
a  (s,(}),  6) 


=   s(x-f)  +  f 

_   r(a+r  cos  <^) 
x-f 


(f ,  (s-l)-sx.) 

1+  -^ , — i_ 


(fg  (S-D-SXq) 


(a+r  cos  ^) 


b(s,ct),0) 
c(s,(}),e) 
d(s,(j),e) 
e(s,(}),e) 


(a+r  cos  (p)   (x-f) 
r 

r(x-f) 


a+r  cos  (}) 


[f  ,  (s-1) -sx  ,  ]  [a+r  cos  (J)] 

<i> $ 


[f „ (s-l) -sx„] r 


a+r  cos  (|) 


The  Euler  equation  for  (4  8)  yields  the  following  equation 
for  g  in  D- 

(49)   1^  (ag^+dg^+eg^)  +  f^Cbg^+dg^)  +  feCcg^+eg^)   =   0 

*    * 

with  boundary  conditions  on  C  ,  T 


(50) 


1^  (a  g^  +  d  g^  +  eg^)   =  0  . 
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4.3.   Difference  Equations  for  the  Vacuum  Region. 

The  usual  approach  for  deriving  difference  equations 
for  g  would  be  to  write  a  discrete  approximation  to  equa- 
tions (49)  and  (50) .   Instead,  we  take  as  our  basic  idea 
the  variational  principle  for  Laplace's  equation,  and  we 
define  a  second  order  accurate  finite  difference  approxi- 
mation to  Dirichlet's  integral  F(g)  given  by  (48).   Then 
F  is  a  function  of  the  values  of  g  at  the  lattice  points, 
and  we  apply  the  variational  principle  requiring  F  to  be 
stationary.   This  last  step  yields  the  difference  equations 
for  g. 

We  define  on  D_  a  lattice  with  uniform  mesh  sizes 
h   ,  h,  ,  hn   and  approximate  the  continuous  function 
g(s,({),e)  by 

^ijk  "  ^^^  ^s'  ^  ^<i>'    ^   ^e^  '   i=0'--"n; 

j=0,...,m;  k=0,...,m 

The  first  step  is  to  write  a  finite   difference 

approximation  to  Dirichlet's  integral.   For  each  rectangle 

in  the  lattice  with  sides   h  ,  h,,  h^  we  make  a  second 

s    (p    9 

order  accurate  difference  approximation  to  the  integrand 
of  (48) ,  and  write  the  integral  as  a  sum  over  all  rectangles 
in  the  lattice.   There  are  many  possible  approximations, 
so  we  choose  one  with  the  minimum  number  of  terms. 

Figure  5  shows  a  rectangle  whose  vertices  numbered 
from  1  to  8  are  points  in  the  lattice.   We  make  the  following 
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second  order  approximation  to  the  integrand  of  (48) 

evaluated  at  the  center  0  of  the  rectangle.   For  any 

fionction  t,  let   t..  ,  i,j  =  1,...,8   denote  the  value 

of   t   evaluated  at  the  middle  point  between  points  i  and  j 

Then 


a  gj  =  I  [(agj)^4  +  Ugl)  ^^   +  (ag^)g7  +  (^^^58^ 
^  5^  =  ?  ^^^21  ^  K^8  -^  (^5^65  -^  (^5')34^ 

2dg^g^  =  f  d^Q[(g^)g^+  ^^s^  58^  f  ^5*^8  +  (5*^56^ 

2eg3g^  =  I  e^^  [(93)32+  (93^7^^(^6)27^  ^59)36] 

+  i^l2f(5sUl  ^  (^3^58^^(^6)45  -^  (^6)18^ 

where  the  derivatives  are  approximated  by  centered 
differences.   We  have  then 

(51)  F(g)  =  h^h^hg  I    (ag^bg2+cg2+2dg^g^+2eg^gQ)Q 

Using  the  previous  approximation  we  have  F  =  F(g.  .,  )  and 
the  variational  principle  implies 

9F 

(52)  T =  0,  1=0,..., n;  j  =  0,...,m;  k=0,...,m 

^^ijk 
which  yields  the  difference  equation  for  each  point  ijk 

in  the  lattice. 
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If   ijk   is  an  interior  point,  (52)  can  be  written  as 

(53)   L,+  Lj+   L^+   L  +  L5=  0,   i=l,...,n-l;  j  =  0,...,m;  k=0,...,ni 
where 

h  =  H;  f(^5s^i+l/2,j,k-  (^^s^i-l/2,j,k^ 
is  a  centered  difference  approximation  to  -5—  (ag  ) 

h  =  t^     f(^5<^^i,j+l/2,k-  (^94>^i,j-l/2,k^ 
is  a  centered  difference  approximation  to  -r-r-  (bg.) 

^3=  K^   f^^^9^i,j,k+i/2-  (^^e^i,j,k-i/2^ 

g 

is  a  centered  difference  approximation  to  r-r  (cg„)  .   The 
derivatives  are  also  approximated  by  centered  differences, 
so  for  example 

^°^e^i,j,k+l/2  "  i^  ^i,j,k+l/2^%,j,k+l"^i,j,k^  • 

L, ,L-,L-   correspond  to  the  usual  difference  approximation 

to  the  corresponding  terms  in  equation  (49) .   The  mixed 

9  9 

derivative  terms   7^—  (dg^)  +  -r-r  (dg  )   are  approximated  by 

dS      (p      d<J)      S 

^4  "  2h~h~  ^'^i+l/2,j+l/2,k^%+l,  j+l,k"%,j,k^ 
~  '^i-l/2,j-l/2,k^^i,j,k~^i-l,j-l,k^ 
^ ^i+1/2 , j-1/2 ,k  ^^i+1 , j-1 ,k"^i jk^ 
*^i-l/2,  j+l/2,k^^i,  j,k"^i-l,  j+l,k^^  ^ 


-36- 


To  interpret  this  approximation  we  make  the  change  of 

variables 

u  =   TT  h    (t—  +   z—) 
2      u  h,         h 

q>  s 

-  =  I  ^<^  -  ^ 

<i         s 
so   that 

(J)    s 

and  we  can  see  that  L.  is  a  centered  difference  approxima- 
tion to  the  right-hand  side.   The  mixed  derivative  terms 

8  9 

g^(egg)  +  9^ (egg)   are  approximated  by 


^5    2h  h„  ^^i+1/2, j,k+l/2^^i+l, j,k+l  ^i, j,k' 


"  ^i-1/2,  j,k-l/2^'^i,  j,k"^i-l,j,k-l^ 

"^^i+1/2, j,k-l/2^^i+l,j,k-l~^i,j,k^ 

^i-1/2, j,k+l/2^^i,j,k"^i-l,j,k+l^^^ 

which  can  be  interpreted  in  a  similar  way.   Equation  (53) 

* 
is  a  15  point  formula  for  Laplace's  equation  in  D-  instead 

of  a  7  point  formula  we  would  get  for  Laplace's  equation 

in  the  original  domain.   This  is  a  consequence  of  the  fact 

that  the  mapping  is  not  con formal  and  therefore  the  equation 

* 
in  D-  contains  mixed  derivatives. 

For  boundary'  points  (52)  can  be  written  as 
(54)   L^  +  I  (^2'^^3^  +  L^  +  L^  =  0,  i=0,...,n;  j,ksO,...,m 
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where  by  definition  the  coefficients  a,d,e  are  zero  outside 
the  domain,  that  is 

^-1/2, j,k  "  '^-l/2,j  +  l/2,k  "  ^-l/2,j,k+l/2  "  ° 
and 

^n+l/2,j,k  "  ^n+1/2, j+l/2,k  "  ®n+l/2 , j ,k+l/2  "  ° 
The  periodicity  conditions  (47)  are  given  by 


^i,m+l,k     ^i,0,k  "^  """1 

"i,-l,k    ~   "i,m,k     1 
(55)  i=0, . . . ,n; 

g.  .   ,,   =   g.  .  ^  +  I~    j,k=0,...,ra 
g.  .   T    =   g .  .    -  I-, 
which  together  with  (53) ,  (54)  determine  the  solution. 


4.4.   variational  Principle   for  Laplace's  Equation  in 
the  Plasma  Region. 

We  follow  the  same  procedure  as  in  the  previous  case. 

* 
We  define  a  mapping  from  D   (plasma  region)  onto  D,  by 


s  = 


f  ((|),e) 
(()  =  <(, 


for  0  <  r  <  f;   0  <  (j)  <  2tt;  0  <  6  <  2Tr;   so  that 
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D,    =    {(s,({),e)     I     Oj<s<_l,    0    <_  (t>    <    2-n;    0£6<2tt} 
r      =    {(s,({),e)     I    s   =    1,    0   £  4)    <_  2tt,    0   £   9    <_  2it} 

Let      4' (r  (s,  <t>,  9)  ,(|),  0)    =   g(s,(|),9).      We  wamt   to    solve 
Laplace's   equation    in   D,     ,    with   Neumann   boundary    conditions 
on    r       ,    and   periodicity   conditions 


(56) 


Vii'di    =    g(s,(i),9+2TT)    -    g(s,<)),e)    =    I, 


Furthermore,  since  g  must  be  harmonic  for  s  =  0 ,  we 
require  that 


(57) 


g(s,<})+2TT,9)  -  g(s,(t),9)  =  0 


The  expression  for  Dirichlet's   integral  has  the  same 

* 

form  as  before  (4  8) ,  where  the  integral  is  taken  over  D, 

and  the  coefficients  are  given  by 


a(s,(}),  9 

b(s,4),9 
c  (s,(t>,  9 

d(s,<^,9 

e  (s,(t),9 


.2  2.2 

f  ,  S    fg 

=   s(a+sf  COS  c}))  (l  +  — 2"  +  2 

f     (a+sf  COS  (^) 

_     a+3f  COS  j) 

S 

sf 


) 


a+sf  COS  (}) 


f  ,  (a+sf  COS  d)) 

_2 ^^ 


2^  4= 

s  f  f , 


a+sf  COS  cf) 


We  define  on   D^   the  same  type  of  lattice  as  before. 

The  equations  for  interior  points  are  the  same  as  before, 

* 

given  by  (53) .   The  free  surface  T        corresponds  to  s  =  1 

and  the  corresponding  'equations  are  given  by  (54)  with  i=n. 
The  periodicity  conditions  (56)  and  (57)  are  given  by 
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i=  0 , .  .  . ,  n ; 
j ,k=0, . . . ,m. 


^i,m+l,k   ^i,0,k 
^i,-l,k   ~  ^i,m,k 


^i,j,m+l    ^i,j,0  "^  ^3 


g .  .   ,   =  g.  .    -  I_ 


4.5.   Laplace's  Equation  at  the  Singular  Curve 

* 

The  points  in   D,   with  s  =  0   correspond  to  points  in 

D,  with  r  =  0,  so  they  are  singular  points  for  the  coordi- 
nate  system.   The  general  solution  in  the  neighborhood  of 
r  =  0  is  of  the  form  k, (6)  +  k_(e)  log  r  +  0(r)   and  since 
we  want  the  solution  to  be  harmonic  at  r  =  0 ,  we  need 
k_(9)  =0.   To  achieve  this,  we  must  formulate  a  mean  value 
theorem  for  k, (9)  in  terms  of  the  values  of  the  solution  at 
adjacent  points. 

If   g   is  harmonic  at   s  =  0,   Laplace's  equation 
implies  that   jr  (bg.)  =  0  at  s  =  0.   Consequently,  we 
eliminate  the  logarithmic  solution  at  s  =  0  by  taking  L2=  0 
(since  L^   is  the  difference  approximation  to  3T(^9(j,))' 

For  any  given  k 

50,j,k  =  g(O,jh^,kh0)  ,      j  =  0,...,m 
is  the  value  of  the  solution  for  r=0,  9=kh.   (i.e.  for 
all  values  of  j   these  represent  the  values  of  the  solution 
at  the  same  point) .   It  then  follows  that  when  we  vary  the 
value  of  the  solution  at  r  =  0 ,  9  =  kh„  ,  we  must  sum  over  j 
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the  variations  due  to  the  variation  of  g^  .  ,  ,  with  the 

0,],k 

condition  that  g-  .is  independent  of  j . 

Therefore,  the  difference  equation  at  a  singular  point 
is  given  by 

I     (L^+  I  L3+L4+L5)^jj^  =  0  ,   i=0,  k=0,...,m 

which  corresponds  to  (54)  summed  over  j,  with  the  additional 
condition  L-  =  0. 

4.6.   Successive  Overrelaxation  Method 


We  solve  the  difference  equations  using  the  successive 

overrelaxation  method  [6].   We  write  the  difference  equations 

in  the  form 

i-1  R 

g^  =  I  \^y^  +  1      ct^g^  "^  ^i  '  i=o,...,R 

k=0        k=i+l 
and  we  define  a  sequence   g.     of  approximations  to  the 
solution  as  follows 

1  1     k=0  ^  ^       k=i+l  ^  ^      ^ 

where   1  <  co  <  2   is  the  relaxation  factor.   Then 

(r\ ) 

g.  =  lim  g.     is  the  desired  solution. 

4.7.   Numerical  Method  for  Updating  the  Free  Boundary  Position. 

The  equation  for  the  free  boundary  T      is  given  by 
r  =  f(4),6),    which  we  approximate  in  the  lattice  by 


:  .-^  =   f(j  h,  ,  k  hg)  ,   j  =  0,...,m;  k=0,...,m 


and  let 
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V=  <^^2  -1^1  -P)jk 

We  denote  by   f .,    the  position  of  the  free  boundary  at  the 

nth   step,  and  by   A  .,    the  value  of   A  .,   evaluated  for 

the  free  boundary   position  at  the  nth  step.   Making  the 

approximation   6v  =  6f  ,  the  path  of  steepest  descent  is 

given  by 

ff"^l^  =   ffj)-  AAf^) 
]k        jk  jk 

This  formula  requires  a  difference  approximation  to  A  .,  , 
which  will  be  discussed  in  the  next  section. 


4.8.   Numerical  Approximation  to  the  Square  of  the  Gradient 
on  the  Free  Boundary. 

Our  experience  shows  that  in  order  to  have  a  stable 
nijmerical  scheme,  it  is  crucial  that  we  choose  an  approxi- 
mation to  the  square  of  the  gradient  on  the  boundary  which 
is  consistent  with  the  discrete  variational  principle  in 
the  interior.   In  principle,  we  could  perform  a  discrete 
interior  variation  of  the  domain,  and  compute  directly  the 

discrete  expression  for  the  first  variation  of  the  energy. 

2 

We  would  get  a  sum  over  the  boundary  which  is  an  0(h  ) 

approximation  to  the  integral  (28) .   This  is  the  proper 
expression  we  should  use  for  our  method  of  steepest  descent. 

This  point  is  made  clear  if  we  compare  the  solutions  for 
a  case  of  an  elliptical  cross  section,  axially  symmetric,  in 
only  two  independent  variables.   In  this  case,  we  know  the 
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solution  to  be  stable,   so  the  scheme  should  yield  solutions 
for  which  the  residual  on  the  boundary  (i.e.  the  maximum 
over  the  boundary  for  the  error  in  the  free  boundary  condi- 
tion (17) )  is  as  small  as  we  want.   Nevertheless,  if  we  use 
an  approximation  to  the  gradient  such  that  the  <t>    and  6 
derivatives  are  approximated  by  centered  differences,  and 
the  s  derivative  by  a  three  point  formula,  we  find  that  the 
residuals  on  the  boundary  will  decrease  to  10~  ,  and  from 
that  point  on  the  scheme  diverges.   We  also  observe  that  if 
we  split  the  mesh  size  in  half,  the  residuals  on  the  boundary 
will  decrease  to  10   ,  and  then  diverge. 

On  the  other  hand,  if  we  use  an  approximation  to  the 

square  of  the  gradient  suggested  by  the  discrete  variational 

-12 
principle,  the  residuals  on  the  boundary  decrease  to  10 

which  is  within  the  round-off  error  for  the  CDC  6600  computer, 

since  in  this  case  the  residual  for  Laplace's  equation  must 

-14 
be  less  than  10 

2 
These  facts  suggest  that  it  is  the  terms  of  0 (h  )  in 

the  approximation  to  the  square  of  the  gradient  on  the 

boundary  which  make  the  numerical  scheme  unstable  in  the 

first  case. 

We  have  chosen  the  expression  below  for  the  square  of 

the  gradient,  as  suggested  by  the  variational  principle. 

We  have  not  carried  out  the  exact  calculation  of  the 

discrete  expression  for  the  first  variation  of  the  energy, 

so  that  the  only  proof  we  have  that  this  approximation  is 
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a  good  one  is  that  it  works  in  practice.   From  (48)  we  have 

that 

2       2      2      2 

(Vtj;)  J  =  ag^  +  bg^  +  cg^  +  2dggg^  +2eg^gg  =  a 

where  J  =  9  (x,y ,  z) /9  (s  ,  <}) ,  9)  .   To  be  consistent  with  the 
variational  principle,  we  approximate  a  in  the  same  manner 
in  which  we  approximated  Dirichlet's  integral,  that  is 

"jk"  4^°'j  +  l/2,k+l/2'^"j  +  l/2,k-l/2'^°'j-l/2,k+l/2"^"j-l/2,k-l/2^ 

and  each  of  the  four  terms  are  approximated  on  the  boundary 
by  the  same  differences  which  were  used  for  the  discrete 
variational  principle.   It  is  important  to  note  that  if  we 
had  evaluated  a.,  directly  instead   of  using  the  previous 
average,  the  resulting  numerical  scheme  would  be  unstable. 


Having   a.,   we  define 
]k 


2     '^ik 
where   J.,  is  the  Jacobian  evaluated  at  the  point  j,k, 
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5.   Results 

We  will  first  discuss  some  of  the  stability  properties 
of  the  numerical  scheme.   It  is  important  to  point  out  that 
the  main  difficulties  arise  from  the  fact  that  we  have  used 
a  formulation  in  terms  of  a  potential  function,  which  leads 
to  Neumann  boundary  conditions.   Since  the  mapping  is  not 
conformal,  in  the  fixed  auxiliary  domain  the  boundary 
conditions  correspond   to  a  linear  combination  of  normal 
and  tangential  derivatives  (50) .   The  same  difficulties 
are  present  if  we  choose  the  same  kind  of  formulation  for 
a  two  dimensional  problem.   Therefore,  we  can  use  the  two 
dimensional  problem  to  check  and  study  the  properties  of 
our  numerical  scheme.   As  it  was  described  in  the  previous 
chapter,  the  main  numerical  problems  were  solved  by  appeal- 
ing to  a  discrete  variational  principle. 

We  have  found  that  the  stability  of  the  successive 
overrelaxation  procedure  for  a  given  position  of  the  free 
boundary  is  strongly  dependent  on  the  ratio  of  the  mesh 
sizes  h  ,h,,h„.   The  method  is  most  stable  and  most 

S   <p   6 

efficient  (i.e.  requires  fewer  iterations)  if  the  mesh 
sizes  are  nearly  equal. 

This  point  is  very  well  illustrated  by  the  following 
comparison.   For  an  axially  symmetric  case  (independent  of  G) , 
we  computed  the  solution  to  Laplace's  equation  with  Neumann 
boundary  conditions  and  fixed  boundaries.   Keeping  constant 
the  total  number  of  points  in  the  mesh,  we  used  mesh  sizes 
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with  h,  =  2TTh   and  h,  =  h  .   In  the  first  case  the  maximum 
cp      s       <p     s 

value  of  the  relaxation  coefficient  03   compatible  with 

convergence  was  1.5  and  it  took  10,000  iterations  to  solve 

-7 
the  equation  within  an  error  of  10   .   In  the  second  case, 

we  could  use  values  of  to  up  to  1.9  and  it  took  only  300 

iterations  to  achieve  the  same  accuracy. 

Another  practical  problem  is  the  computation  of  the 

inductance  matrix.   For  a  given  position  of  the  boundary, 

we  have  to  allow  for  a  given  error  in  the  solution  of 

Laplace's    equation    in    the    interior,    which   in    turn  will 

produce  an  error  in  the  inductance  matrix.   However,  we 

can  minimize  this  error  if  we  use  (32)  to  compute  the 

matrix,  instead  of  (30) .   The  reason  is  that  the  volume 

integral  given  by  (32)  is  stationary  for  solutions  of 

Laplace's    equation,    so    that   the   error   in    this    case  will   be 

of  second  order.   Of  course  this  will  require  an  extra 

amount  of  computation,  but  it  is  worthwhile  since  in  practice 

we  will  perform  several  iterations  in  the.  interior  (of  the 

order  of  four  to  ten)  for  each  iteration  of  the  boundary. 

5.1.   Description  of  the  Physical  and  Numerical  Parameters. 

The  physical  parameters  to  be  prescribed  for  a  given 
geometry  are  the  three  fluxes  and  the  mass  of  the  plasma. 

Instead  of  prescribing  the  mass  of  the  plasma,  we 
prescribed  the  initial  position  of  the  free  boundary  (thus 
prescribing  the  initial  volume  of  the  plasma)  and  the  initial 
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value  of  the  pressure  in  the  plasma.   From  these  two  values 
we  computed  the  mass  of  the  plasma,  using  an  equation  of 
state  with  y   =   5/3. 

We  also  chose  to  prescribe  the  initial  value  of  the 
currents   IwIp^I^  ,    and  we  computed  the  corresponding 
fluxes   FwF„,F^   for  the  initial  position  of  the  free 
boundary.   Once  the  fluxes  and  mass  had  been  computed  for 
the  initial  position  of  the  free  boundary,  they  were  kept 
fixed  throughout  the  rest  of  the  computation. 

As  one  of  the  results  we  have  computed 


6  = 


2^2 


r 

the  average  over  the  free  boundary  of  the  ratio  of  the 
plasma  pressure  to  the  magnetic  pressure  in  the  vacuum. 
A  better  parameter  which  we  can  compute  would  be  the  ratio 
of  the  internal  energy  E^  to  the  total  energy  E  of  the 
system,  but  we  will  use   6   since  it  is  widely  used  in 
the  literature . 

We  will  define  two  important  numerical  parameters. 
First,  rl   2    12 


e„   =  max 


2  ^2  -  2  ^1  -  P 


T   -  "7"        1  „2 
^         2  ^2 

the  absolute  value  of  the  maximum  error  in  the  equilibriiom 
condition,  and  second, 

e^  =    max   max  |  Aip  .  | 

°    i=l,2,3   D. 
D 

the  maximum  error  for  Laplace's  equation  in  the  interior. 
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where  D.  is  the  corresponding  domain  for  \p .  .      The  minimax 

characteristic  of  the  problem  requires  that  we  solve 

Laplace's  equation  very  accurately  in  order  to  satisfy 

the  free  boundary  condition.   On  the   other  hand,  it  is 

a  waste  of  computation  to  make   e^  very  small  when  we  are 

still  far  away  from  the  equilibrium  position.   In  practice, 

we  required   e    to  be  small  compared  to  e^  ,  so  that  we 

prescribed  -:  .  , 

-4 
e   =  min  (10   ,  0.1  £p/c,) 

where   min(a,b)   denotes  the  minimum  between  a  and  b. 
The  choice  of  the  scale  factor  c^  is  due  to  the  fact  that 
the  biggest  error  in  Laplace's  equation  occurs  in  the 
potential  ii,    ,   which  must  be  multiplied  by  the  factor  c,  . 

The  various  convergence  factors  which  occur  in  the 
formulation  of  the  problem  were  determined  by  trial  and 
error  to  produce  a  convergent  scheme.   For  the  relaxation 
factor  to   in  the  successive  overrelaxation  method  we  used 
1.9  for  the  vacuum  region,  and  1.8  for  the  plasma  region. 
For  the  relaxation  factor  X    in  the  steepest  descent  method 
we  used  values  ranging  from  0.02  to  0.06. 

5.2.   Description  of  the  Mesh  Size  and  Computing  Time. 

For  our  production  runs  we  have  used  a  mesh  with 
30  points  in  the  <}>  and  9   directions  and  12  points  in  the 
radial  direction,  six  of  which  lie  in  the  vacuum  region 
and  the  other  six  in  the  plasma  region.   The  program  was 
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designed  by  trying  to  minimize  the  core  space  needed.  For 
the  mesh  size  just  described,  the  amount  of  core  needed 
is  110  K  which  is  one-third  of  the  total  core  available 
in  the  CDC  6600  computer.   A  program  which  takes  less  time 
per  iteration  can  be  written  using  a  larger  core  space, 
but  for  practical  purposes  this  is  not  convenient. 

The  computing  time  for  a  solution  in  this  kind  of  a 
mesh  varies  with  the  geometry  and  stability  of  the  solution. 
For  a  stable  configuration  it  takes  around  40  minutes  to 
yield  a  solution  with   £„  =  10   .   In  this  case,  1800  itera- 
tions were  needed  to  solve  Laplace's  equation  for  300 
iterations  of  the  boundary  position. 

We  also  ran  some  cases  in  a  bigger  mesh  to  study  the 
dependence  of  the  solution   in  the  mesh  size.   We  started 
with  a  mesh  with  25  points  in  the  (j)  and  9  directions,  and 
10  points  in  the  radial  directions.   Then,  we  split  the 
mesh  size  in  half,  getting  a  mesh  with  50  points  in  the 
(|)  and  9   directions,  and  18  points  in  the  radial  directions, 
This  gave  us  2  500  points  in  the  free  boundary,  but  it 
required  330  K  of  core  space,  which  is  the  total  capacity 
of  the  CDC  66  00.   In  this  case,  three  hours  of  computing 
time  were  required  to  yield  a  solution  with  Cp  =  10    for 
a  stable  case. 

5.3.   Axially  Symmetric  Torus  with  Circular  Cross-Section. 

As  a  first  test  of  our  model,  we  have  studied  the 
equilibrium  and  stability  properties  of  an  axially 
symm.etric  torus  with  circular  cross-section.  Since  we  are 
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specifically  interested  in  the  Tokamak,  we  have  used  physical 
parameters  comparable  to  those  of  the  Tokamak.   We  are 
mainly  interested  in  studying  the  stability  properties 
as  we  vary  the  aspect  ratio;  i.e.,  the  ratio  of  the  major 
radius  to  the  minor  radius  of  the  torus. 

Our  results  are  similar  to  those  given  by  Richtmyer 
et  al.  in  [9],  that  is,  a  torus  with  a  small  aspect  ratio 
is  more  stable  than  one  with  a  larger  aspect  ratio. 

We  have  studied  in  particular  the  cases  with  ratios 
10:4  and  10:3.   In  practice,  for  technical  reasons  the 
lowest  aspect  ratio  that  can  be  achieved  in  a  Tokamak 
is  around  3,  so  that  the  second  case  corresponds  to  a 
ratio  which  is  found  in  practice. 

For  a  fixed  geometry,  we  studied  the  dependence  of 
stability  properties  on  the  physical  parameters.  As 
expected,  we  found  that  an  increase  of  the  toroidal  field 
(i.e.,  Bq)  will  make  the  configuration  more  stable.  On  the 
other  hand,  if  we  keep  the  toroidal  field  fixed  and  increase 
the  plasma  current  (i.e.  the  poloidal  field),  at  some  point 
the  configuration  will  become  unstable.   Furthermore,  if  the 
plasma  current  tends  to  zero,  no  equilibrium  exists  and 
the  plasma  goes  all  the  way  to  the  outer  wall. 

We  measure  the  stability  of  a  configuration  by  looking 
at  the  convergence  of  the  free  boiondary  condition.  We  can 
catalog  the  results  in  three  different  groups.   The  first 
one  is  the  case  for  which  £„  converges  to  zero  which  indi- 
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cates  that  the  energy  has  a  minimum  at  the   equilibrium 
position,  and  therefore  it  corresponds  to  stable  equilibrium. 
The  second  case  for  which  e„   diverges      corresponds 
to  unstable  equilibrium  or  to  the  fact  that  equilibrium 
just  doesn't  exist.   In  this  instance  the  energy  decreases 
for  each  boundary  position  until  the  plasma  touches  the 
outer  wall.   In  between  we  have  a  whole  variety  of  cases 
for  which   £„   stays  bounded  but  it  does  not  converge  to 
zero;  and  the  free  boundary  oscilates  around  the  equilibrium 
position  with  amplitudes  roughly  proportional  to  e„. 

We  believe  that  the  explanation  for  this  behavior   lies 
with  the  nonlinearity   of  the  free  boundary  condition.  Due 
to  a  bifurcation  process  more  than  one  equilibrium  position 
exists,  and  as  the  axially  symmetric  equilibrium  becomes 
unstable  the  numerical  method  is  not  accurate  enough  to 
distinguish  between  neighboring  solutions. 

In  Table  I,  results  for  an  aspect  ratio  of  10:4  are 
shown.   The  initial  position   of  the  free  boundary  is  a 
torus  of  minor  radius  equal  to  3.   The  initial  value  of  the 
currents  are  given  by  1^,1^,1^        as  defined  in  Chapter  2. 
Physically  I,  corresponds  to  the  plasma  current.  The 
pressure  is  denoted  by  p,  and  6,  Sp   were  defined  in 

Section  5.1.   For  stable  cases,  the  error  e^  decreases 

-7 
rapidly  to  10   ,  at  which  point  we  stopped  the  process. 

-7 
Therefore  when  reading  the  results  ep  =  10    corresponds 

to  a  stable  case. 
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In  the  first  four  cases  denoted  by  N  =  1,2,3,4,  we  have 
kept  fixed  the  poloidal  field,  the  plasma  pressure  and  we 

varied  the  toroidal  field.   We  see  that  the  first  two  cases 

-7 
are  stable  as  shown  by  e„  =  10   .   In  cases  3  and  4,  a  very 

mild  instability  appears  as  shown  by  the  value  e^  =  10 

This  means  that  the  process  will  not  converge   any  further 

—  6 
than  10   ,  but  on  the  other   hand  it  will  not  diverge  and 

the  free  boundary  oscillated  around  the  equilibrium  position 

with  extremely  small  amplitudes. 

In  cases  5  to  8,  we  have  fixed  the  values  of  the 
currents,  but  we  have  changed  the  values  of  the  initial 
pressure  (i.e.  the  mass  of  the  plasma).   Cases  5  to  7  are 
stable,  and  case  8  has  a  very  mild  instability. 

In  general  we  find  that  for  this  aspect  ratio, 
configurations  with  8  up  to  0.05  are  stable,  and  configura- 
tions with  3  up  to  0.09  are  very  nearly  stable. 

In  case  number  9  the  plasma  current  is  extremely  small, 
and  in  this  case  an  equilibrium  configuration  does  not  exist 
and  the  plasma  expands  to  the  outer  wall. 

Figure  5.1  shows  the  equilibrium  position  for  case  1 
of  Table  I.   In  this  case  the  equilibrium  configuration  is 
axially  symmetric,  and  its  cross-section  is  a  circle  with 
its  center  displaced  away  from  the  axis  of  the  torus. 

In  Table  II,  results  for  an  aspect  ratio  of  10:3  are 
shown.   The  initial  position  of  the  free  boundary  is  a  torus 
of  minor  radius  equal  to  2 . 
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In  cases  1  to  7   we  kept  fixed  the  plasma  pressure 
and  the  plasma  current,  while  we  decrease  the  toroidal 
field.   We  can  see  that   £„   increases  as  the  toroidal 
field  decreases,  so  that  the  configuration  becomes  more 
unstable.  In  cases  5,  9  to  14   we  have  kept  fixed  the 
toroidal  field  and  we  have  increased  the  plasma  current 

and  the  plasma  pressure   in  such  a  way  that  the  initial 

1   2 
pressure  is  the  average  of  -j   B,   (poloidal  field)  over 

the  free  boundary.   These  results  indicate  that  the  con- 
figuration  becomes  more  unstable  as  we  increase  the 
plasma  current.   Cases  4,  15  to  17   show  similar  results 
for  a  different  value  of  the  toroidal  field. 

In  general  we  find  that  even  for  values  of  B  as  small 
as  0.02   the  configurations  are  not  stable.   If  we  compare 
the  results  in  Table  I  and  Table  II,  we  conclude  that  for 
similar  values  of  6   a  torus  with  aspect  ratio  10:4  is  more 
stable  than  a  torus  with  aspect  ratio  10:3. 

5.4.  Axially  Symmetric  Torus  with  Elliptical  Cross-Section. 

One  important  question  which  has  been  raised  lately  is 
that  of  the  stability  properties  of  an  axially  symmetric 
torus  with  a  noncircular   cross-section.   In  particular  we 
want  to  compare  the  stability  properties  of  an  elliptical 
cross-section  with  those  of  a  circular  cross-section. 

In  Table  III,  results  for  an  elliptical  cross-section 
are  shown.   The  ratio  of  the  major  axis  to  the  minor  axis  a 
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of  the  ellipse  is  1.5.  Figure  5.2  shows  the  equilibrium 
position  for  case  2  in  Table  III.   In  this  case  the  equili- 
brium configuration  is  axially  symmetric,  and  the  cross- 
section  is  an  ellipse  whose  center  has  been  displaced  away 
from  the  axis  of  the  torus. 

From  Table  III  we  can  draw  the  following  conclusions. 
Case  1  is  stable,  as  we  would  expect  from  its  10:4  aspect 
ratio.   Cases  2  to  4  show  that  the  configuration  becomes 
more  unstable  as  we  decrease  the  toroidal  field.  A  comparison 
of  cases  2  and  3  with  cases  6  and  7  shows  that  the  stability 
becomes  much  worse  as  we  increase  the  aspect  ratio.  Note 
that  cases  2  and  6  correspond  to  the  same  value  of  6,  but 
a  small  increase  in  the  aspect  ratio  (from  10:3  to  10:2.45) 
makes  a  big  difference  in  the  value  of  e^  . 

The  conclusion  is  that  the  dependence  of  stability 
properties  on  the  physical  parameters  and  the  aspect  ratio 
is  similar  for  an  elliptical  cross  section  to  that  of  a 
circular  cross  section. 

Using  Tables  II  and  III  we  can  compare  the  stability 
of  an  elliptical  cross  section  with  that  of  a  circular  cross 
section,  for  similar  physical  parameters  and  aspect  ratio. 
We  see  that  both  yield  similar  results,  and  that  the 
stability  properties  are  comparable  for  equal  values  of  3. 
We  conclude  then  that  our  model  shows  that  an  elliptical 
cross  section  and  a  circular  cross  section  have  essentially 
the  same  stability  properties. 
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5.5.   Configurations  without  Axial  Symmetry. 

We  have  computed  equilibrium  configurations  for  several 
types  of  three  dimensional  geometries  in  order  to  compare 
their  stability  properties  with  those  of  axially  symmetric 
cases. 

The  results  shown  in  Table  IV  correspond  to  a  torus 
with  circular  cross  section  whose  center  describes  a  helix 
as  9  varies  in  2Tr. 

In  cases  1  to  9  we  have  kept  the  physical  parameters 
constant,  while  we  varied  the  geometry.   Case  1  corresponds 
to  the  axially  symmetric  torus   and  by  comparison  with  the 
other  cases  we  see  that  the  effect  of  introducing  a  helical 
perturbation  in  the  outer  boundary  is  to  make  the  configura- 
tion more  stable.   The  same  conclusion  holds  true  if  we 
compare  cases  9  and  10.   Case  11  is  stable  as  we  would  expect 
from  its  10:4  aspect  ratio. 

Figures  5.3  and  5.4  show  the  outer  boundary  and  free 
boundary  for  case  11  in  Table  IV.  We  have  plotted  the  surfaces 
r  =    f{<\)  ,d)      in  the  coordinate  system  defined  in  Chapter  4. 
It  is  useful  to  remember  that  a  plane  with  r  =  constant 
corresponds  to  an  axially  symmetric  torus  with  circular 
cross  section.   From  these  figures  we  observe  that  the 
solution  for  the  free  boundary  is  smooth  and  that  it  takes 
the  same  kind  of  shape  as  the  outer  boundary.  By  looking  at 
the  dependence  on   4)  we  can  see  that  the  plasma  surface 
has  been  displaced  away  from  the  axis  of  the  torus. 
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Figures  5.5  and  5.6  show  the  outer  boundary  and  free 
boundary  for  case  4  in  Table  IV   We  can  clearly  see  the 
four  helical  turns  of  the  outer  surface  and  the  free 
boundary. 

The  results  shown  in  Table  V  correspond  to  a  perturba- 
tion of  the  type  A  cos  (k-fJi+k-S).   The  results  show  that  for 
6  =  0.09,  cases  3,4,5  are  more  stable  than  the  axial ly 
symmetric  case,  and  for  6  =  0.2  8  there  is  no  improvement. 

Figures  5.7  and  5.8  show  the  outer  boundary  and  free 
boundary  for  case  4  in  T^ble  V,  and  figures  5.9  and  5.10 
show  the  outer  boundary  and  free  boundary  for  case  5 
in  Table  V. 

The  results  shown  in  Table  VI  correspond  to  an  outer 
surface  which  is  obtained  by  rotating  an  elliptical  cross 
section  around  its  center  as  0  varies  in  2tt.   They  show 
that  the  three  dimensional  configurations  obtained  are 
more  stable  than  the  axially  symmetric  case. 

5.6.   Dependence  of  the  Solutions  on  the  Mesh  Size. 

To  study  the  dependence  of  the  solution  as  we  vary 
the  mesh  size,  we  ran  some  of  the  previous  cases  for  the 
biggest  possible  mesh  within  the  capacity  of  the  CDC  6600 
computer  (described  in  Section  5.2).   We  found  that  changing 
the  mesh  size  did  not  change  the  kind  of  solutions  obtained. 
That  is,  those  cases  which  were  staJale  for  a  smaller  mesh 
were  also  stable  for  the  bigger  mesh  and  those  cases  which 
were  imstable  yielded  the  same  values  of  Cp  as  previously. 
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2 
We  expect  the  discrete  solution  to  converge  as  0 (h  ) 

to  a  continuous  solution,  as  h  tends  to  zero,  since  all 

2 
the  approximations  made  were  0(h).   To  check  this  it  would 

be  necessary  to  split  the  mesh  size  by  half  at  least  two 

times,  which  we  cannot  do  because  of  the  limited  capacity 

of  the  computer.   The  existence  of  a  continuous  solution 

to  the  equilibrium  problem  is  another  matter,  and  all  we 

can  say  is  that  since  the  discrete  solutions  are  quite 

smooth  we  would  expect  them  to  converge  to  a  continuous 

solution  as  h  tends  to  zero. 

5.7.   Conclusions. 

The  main  conclusion  is  that  we  have  developed  a 
numerical  scheme  to  compute  configurations  of  fully  three 
dimensional  magnetohydrodynamic  equilibrium,  which  also 
allows  us  to  study  their  stability  properties. 

Our  results  show  that  a  torus  with  small  aspect  ratio 
is  more  stable  than  one  with  a  large  aspect  ratio.  For  the 
same  aspect  ratio,  an  axially  symmetric  torus  with  circular 
cross  section  and  an  axially  symmetric  torus  with  elliptical 
cross  section  have  similar   stability  properties. 

We  also  found  several  three  dimensional  configurations 
which  are  more  stable  than  neighboring  axially  symmetric 
cases.   However,  the  improvement  in  the  stability  properties 
by  introducing  a  three  dimensional  perturbation  of  the  outer 
boundary  is  not  nearly  as  big  as  the  improvement  obtained 
by  lowering  the  aspect  ratio.  This  might  not  be  the  case 
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for  three  dimensional  configurations  which  are  very  different 
from  the  axially  symmetric  case,  a  possibility  which  we  can 
explore  by  extending  our  method  to  a  more  general  coordinate 
system  as  previously  discussed. 
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TABLE  I 
Axially  Symmetric  -  Circular  Cross  Section 
Big  Radius  10    Small  Radius  4    Plasma  Radius  3 


No. 

h 

^2 

^3 

2p 

6 

Cp 

1 

10 

180 

180 

0.30 

0.030 

10-^ 

2 

10 

150 

150 

0.30 

0.041 

10-^ 

3 

10 

130 

130 

0.30 

0.054 

10-^ 

4 

10 

110 

110 

0.30 

0.073 

10-^ 

5 

10 

180 

180 

0.10 

0.010 

10-^ 

6 

10 

180 

180 

0.25 

0.024 

10-^ 

7 

10 

180 

180 

0.36 

0.034 

10-^ 

8 

10 

180 

180 

1.00 

0.086 

10-^ 

9 

1 

180 

180 

0.36 

0.030 

00 
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TABLE  II 
Axially  Symmetric  -  Circular  Cross  Section 
Major  Radius  10    Minor  Radius  3    Plasma  Radius  2 


No. 

h 



^2 

^3 

2p 

6 

ep 

1 

10 

270 

270 

0.65 

0.03 

10-^ 

2 

10 

230 

230 

0.65 

0.04 

10-^ 

3 

10 

200 

200 

0.65 

0.05 

10-^ 

4 

10 

180 

180 

0.65 

0.06 

10-^ 

5 

10 

150 

150 

0.65 

0.09 

10-2 

6 

10 

130 

130 

0.65 

0.12 

1 

7 

10 

80 

80 

0.65 

0.27 

00 

8 

10 

149 

151 

0.65 

0.09 

10-2 

9 

10 

151 

149 

0.65 

0.09 

10-2 

10 

15 

150 

150 

1.47 

0.18 

10-^ 

11 

20 

150 

150 

2.60 

0.28 

1 

12 

25 

150 

150 

4.08 

0.38 

1 

13 

30 

150 

150 

5.87 

0.47 

1 

14 

40 

150 

150 

10.40 

0.60 

oo 

15 

6 

180 

180 

0.23 

0.02 

10-^ 

16 

8 

180 

180 

0.41 

0.04 

10-^ 

17 

15 

180 

180 

1.46 

0.13 

1 
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TABLE    III 


Axially   Synunetric    -   Elliptical   Cross    Section   - 


1.5 


Major  Radius  10 


No. 

Outer 
Shell 
a 

Plasma 
a 

^1 

^2 

^3 

2p 

B 

Ep 

1 

4 

3 

10 

180 

180 

0.114 

0.01 

10-^ 

2 

3 

2 

15 

180 

180 

0.554 

0.05 

10-^ 

3 

3 

2 

15 

150 

150 

0.554 

0.08 

10-^ 

4 

3 

2 

15 

120 

120 

0.554 

0.12 

1 

5 

3 

2 

20 

120 

120 

0.980 

0.19 

1 

6 

2.45 

1.63 

10 

150 

150 

0.365 

0.05 

10-^ 

7 

2.45 

1.63 

10 

130 

130 

0.365 

0.07 

1 
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TABLE  IV 

Major  Radius  10    Minor  Radius  RO    Plasma  Radius  RI 

2        2  2  ? 

r  ((J),6)=R0  [  (cos(j)+Acos  k-Scos  k-S)  +  (sin({)+Acos  k,9sin  k20)  ] 


No. 

RO 

RI 

A 

^ 

^2 

h 

^2 

^3 

2p 

6 

^r 

1 

3 

2 

0.0 

- 

- 

10 

150 

150 

0.65 

0.09 

10-2 

2 

3 

2 

0.10 

0 

2 

10 

150 

150 

0.65 

0.09 

10-^ 

3 

3 

2 

0.10 

0 

2 

10 

150 

150 

0.65 

0.09 

10-^ 

4 

3 

2 

0.15 

0 

4 

10 

150 

150 

0.65 

0.09 

10-^ 

5 

3 

2 

0.15 

1 

1 

10 

150 

150 

0.65 

0.09 

10-^ 

6 

3 

2 

0.15 

2 

2 

10 

150 

150 

0.65 

0.09 

10-^ 

7 

3 

2 

0.15 

4 

4 

10 

150 

150 

0.65 

0.09 

10-^ 

8 

3 

2 

0.15 

6 

6 

10 

150 

150 

0.65 

0.09 

10-^ 

9 

3 

2 

0.0 

- 

- 

20 

150 

150 

2.61 

0.28 

1 

10 

3 

2 

0.15 

0 

4 

20 

150 

150 

2.61 

0.28 

10-1 

11 

4 

3 

0.10 

1 

1 

10 

180 

180 

0.30 

0.03 

10-^ 

12 

2 

1.33 

0.10 

2 

2 

7 

180 

180 

0.71 

0.07 

10-3 

13 

3 

2 

0.15 

2 

2 

10 

180 

180 

0.65 

0.06 

10-^ 

-64- 


TABLE  V 


Major  Radius  10    Minor  Radius  3    Plasma  Radius  2 
r((|),e)  =  3(1.0  +  A  cos  (k^4)+k2e)) 


No. 

A 

''l 

^2 

h 

^2 

^3 

2p 

6 

-r 

1 

0.0 

- 

- 

10 

150 

150 

0.65 

0.09 

10-2 

2 

0.1 

0 

1 

10 

150 

150 

0.65 

0.09 

1 

3 

0.1 

0 

2 

10 

150 

150 

0.65 

0.09 

10-^ 

4 

0.1 

1 

2 

10 

150 

150 

0.65 

0.09 

10-^ 

5 

0.1 

1 

4 

10 

150 

150 

0.65 

0.09 

10-^ 

6 

0.0 

- 

- 

20 

150 

150 

2.61 

0.28 

1 

7 

0.1 

0 

8 

20 

150 

150 

2.61 

0.28 

1 

8 

0.1 

1 

4 

20 

150 

150 

2.61 

0.28 

1 

9 

0.1 

2 

4 

20 

150 

150 

2.61 

0.28 

1 
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TABLE  VI 

The  outer  surface  is  obtained  by  rotating  the  elliptical 
cross  section  as  we  go  around  the  torus.  The  cross  sec- 
tion is  rotated   k    full  turns  as    6    varies  in  2tt. 


Major  Radius  10  Elliptical  Cross  Section  —  =  1.5,  a  =  3 

3. 


No. 

k 

h 

^2 

^3 

2p 

e 

^r 

1 

0 

15 

120 

120 

0.55 

0.12 

1 

2 

1 

15 

120 

120 

0.55 

0.12 

10-1 

3 

2 

15 

120 

120 

0.55 

0.12 

10-1 

4 

4 

15 

120 

120 

0.55 

0.12 

10-2 
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fig.  I 
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fig.  2 
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fig.  3 
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fig.4 
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HETfl   =  324.0 


Fig.  5.1.  Axially  Symmetric  Torus  with  Circular  Cross  Section, 
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HEin   =  324.0 


Fig.  5.2.  Axially  Symmetric  Torus  with  Elliptical  Cross  Section, 
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c 

c 


PROGRAM  BETA( I MPUT, OUTPUT, TAPEl , TAFE2 ) 

COMMON  G(9,5  2,52),A{9,52.52),U(9',52,52),R(52;5?),Fi5?,«2)|T(52i52) 
liCl(9,52)iC2(9,52),El(9,52),E2(9',5Z)'.Cr(52),Pl(53):B2(«2>,Dl(52), 
2D2(5I2),RA,RC,RI,NI  ,  N  J,  PETE  •  PEF  I  ,  BE^  E  .  GRA  i  H  ,  HS  ,  VOL  .  ENER',  RELB 

COMMON  EX(52,52),EF(52.52)iET(52'.52).A11,A12:a22,A33 

lALl,«LOCr(A33)-LOCr(G(l,l,l))*l 

CALL  ZER0(G, IALL) 

READ  IN  INITIAL  PARAMETERS 

RfcAO  1,RA,RC,RI 

1  F0RMAT(3F7,2) 

RfcAO  2,PETE,PEFI,BETE,QRa 

2  rORMAT(4r7,2) 
R£AO    3,NI«NJ 

3  rORMAT(2M) 

REAU  4,RELI,RULA,RELB,ERRIiERRB 

4  F0RMAT(3F6.3,2F10,7) 

PRINT  5,RA,RP,RI 
5   F0RMAT(><l^,2y,^BIG  RAD  lUS^i  2X,F7  .  2i  2X,  ^SMALL  RAD  lUS^.  2y  ,  F7  ,  2,  2X, 
IJ'INITIAL  SURFACE^, 2X.F7,?) 
PRINT  6,PETE,PFFI,BETE.GRA 

6  FQRMAT(^  J«i2Xi)'VACUUM  THFTA  period', 2X,r7, 2, 2X,^VASUUM  Fl  PERIOD?!, 
12X,r7,  2,  2X,  ^PLASMA  THETA  PER  I  OD^',  2)< ,  F7  ,2.  2X,  *  I  M  T  I  AL  PPESSUREj<i  2X, 
2F7,2) 

PRINT  7,RELI.RULA,RELB,ERRI,ERRB 

7  FORMAT(^  ^,2X, ^VACUUM  REL^ i 2X , F6 . 3 1 2X, ^PLASMA  REL^T2X, F6  ,  3, 2X, 
l?!bOUND  REl^,2X,F6, 3, 2X. ^INTERIOR  ERRORj<,  2X,  Fj  ,  3  ,  2X7^001  NO  ERROR 
2?!i2X,riO,7) 

R6A0  1557,ASyE,FAC 

1557  F0RMAT(2F7,5) 

REAU  1558,FEC,0RRI,0RRB 

1558  F0RMAT(F5,3,2Fl5.7) 

PRINT  1559,FaC,FEC,0RRI,0RRB 

1559  FORMAT(;<  ^  ,  3X,  2F10  ,  3,  2F15  ,  7  ) 
GAS|5, 0/3.0 

READ  1580# INn,RCAT,TOM 

1580  FOR«AT( l2,Fln,7,F10,l) 

IF  IND  LESS  THAN  0,WE  USE  AS  INPUT  PREVIOUS  S0LUTI8N  STORED 

IN  TAPEl, 

1F(  IND.LE.O)  GO  TO  1981 

COMPUTE  AXIALLY  SYMMETRIC  SOLUTION 

CALL  ASYM(FLU1,FLU2,FLU3,ASYE) 

C0Ns0,5*GRA*(VriL**6AS) 

1581  CONTINUE 
UlNiO.O 
ARRJsERRi 
BMAX'1.0 
Pl=iil4l5926536 
PP  =  1,0/(2,0*PI  ) 
H"(2,0*PI )/(Nj»0.0) 
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1560 


1561 
2000 


710 


INITIAL  GUESS  FOR  FREE  BOUNTARY 


HS=1,0/(NI-1.0) 

Nl«NJ*l 

N25NJ*2 

N3k2 

R6r§l,0 

HA»0,5/HS 

HBbO,5/H 

HC«(H*H)/(HS*HS) 

HO«H/(2,0*HS) 

Rfc08l,0-RELI 

RAOfl.O-RULA 

NllBNI-1 
RA2iO,5*RA 

T6R81, 0/3,0 

IF"(  JND.LT.O)  Gn  TO  1560 

COMPUTE  OUTER  BOUNDARY  AND 

CALl.  SURF 

CALL  SARF 

S6T  INITIAL  SOLUTION  EQUAL  TO  SCLuTlON  TO  AXfALLV  SYMMFTRIC  CASE 

DO  a  J«1,N2 

DO  d  K«1,N2 

DO  8  I«1,NI 

A( I| J,K)«PP*(K-2)*W 

U(  I|  JiK)tiPP*(K-2)*M 

G{1| J,K)«C1( !# J) 

GO  fO  1561 

REWJND  1 

Njs4 

READ  INITIAL  SOLUTION  FROM  TAPEl 

REAUd)  (<  (A(  I,  J,K),Jsl,N2)#Ksl,N2), 

( < (L( I, J.K), J=i,N2>*Ksi,N2) 

(((G(I, J,K), j8i,N2)#Ksi,N2) 

((R(J.K), J.l,N2).K»l,N2) 

( (EX( J.K), Jal,N2)iKsl,N2) 

FLU1,FLU2,FLU3,C0N 


laliNI) 
,ial,NI) 
:l,NI) 


l! 


READ(l) 

REAOd) 

REA0<1) 

R6AU<1) 

REAOd) 

CONTINUE 

CONTINUE 

DO  9  J«1,N2 

CF( J)«COS( ( J-2)*H) 

COMPUTE  DERIVATIVES  OF  OUTER  SURFACE 

DO  no  K«2,M 

DO  710  jn2,M 

EF( JiK)sHB*(EX( J*1,K)«EX(J"1»K)  > 

ET<JiK)sHB*(EX(J,K*l).EX(JiK-l)) 

PkRIODlClTY  CONDITIONS  FOR  DERIVATIVES  OF  OUTER  SURFACP 

DO  711  K»1,N? 

EF(l,K)aEF(M,K) 

ET{liK)sET(M,K) 
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Er(N2,K)aEr(?,K) 

711  ET{N2,K)sET(?,K) 
DO  712  J«liN2 
EF{J,1>»EF(J,N1) 
ET(J,l)sET(j,Nl) 
EF(J,N2>»Er(J,2) 

712  ET(J,N2)«ET( J,2) 
ITE«»0 

NTEH-O 
INEH»0 
PRINT  10 

10  rORMAT(^l><,4X,)<ITERATION  NUMBER^',  5^  ,?<  VOLUME?!;  5)<,  ^PRESSl  RE^,  8X, 
li<MAX  RES^i7X,^MAX  RES  I  NTER I  Or^  ,  7X  i  ^MAX  ERROR  CN  B9UNDO 

PERIODICITY  CONDITIONS  FOR  FREE  BOUNDARY 
100   DO  603  K«1,N? 

R<1|K)bR(N1,K) 

R{N2iK)sR(2,K) 
604   CONTINUE 

DO  604  J«1,N2 

R(J|1)bR(J,M) 

R(J|N2)3R(J,2) 
604   CONTINUE 

COMPUTE  DERIVATIVES  OF  FREE  BOUNDARY 

DO  11  Ks2|Nl 

DO  11  J82|N1 

F<J|K)bHB*(R( J*1,K).R(J-1,K)) 

T( J|K)bHB*(R( J,K*1)-R(J,K«1)) 

11  CONTINUE 

PERIODICITY  CONDITIONS  FOR  DERIVATIVES  OF  FREE  BPUNDARV 

DO  600  K-liN2 

F<1|K)«F(N1,K) 

T(1|K).T(N1|K) 

F(N2iK)bF(2,K> 

T(N2.K)aT<2iK) 

600  CONTINUE 

DO  601  J«1,K2 

F(Jil)«r(J#M) 

T(Jjl>«T{J.M) 

F(J,N2).F(J,?) 

T(J|N2)bT(J,2) 

601  CONTINUE 

ITERATE  SOLUTION  FOR  POTENTIALS  G  And  A  IN  Ti-E  VACUUM  PEGION 
200  UbKO.O 

DO  72  lsl,Nl 

SV«{  N0,5)*hS«'1.0 

DO  71  J82#N1 

Rl«0i5*(R(J»2)*R(J.Nl)) 

R280,5*(T(J,2>*T(J,N1)) 


-83- 


RiJBQ,5*(EX(j,2)*EX(J.Nl)  ) 

R4b0.5*(ET(J,2)*ET(J,N1)) 

R3sn"1.0)*hS*(R5-Rl)*Ri 

Cl(  li  JjsCRS-CR'S.Rl)  )/(RA*R3*Cr(J)  ) 

R38{I.0,5)*HS*(R5«R1)*R1 

El(|iJ)«(HD*R3*(SV*R2-R4*(  I"0,5)*HS)  )/(RA*R3*CF(J)^ 
7i   CONTINUE 
72   CONTINUE 

DO    i2    K=2iNl 

R4=U,5*(Cr(2)*CF(Nl)) 

DO    70    I«1,NI 

Rl3U,5*(R(2,K)*R(Nl,K)) 

R2sQ,5*(F(2jK)*F(Nl.K)) 

R6aU,5*(EX(2,K)*FX(Nl,K)) 

RbsOi5*{EF(2,K)*EF(Nl,K>) 

Sf(|-1)*HS 

R38S*(R6-Rl)*Rl 

BK I )»(RA*R3*R4)*(R6-R1)/R3 

R3s(S*0,5*HS)*(R6^Rl)*Rl 

DX( pa(HD*(RA*R3*R4)*(R2*(S*0,5*HS"1.0>-(S*0l5*HS)»R5))/R3 
70      CONTINUE 

DO    13    Jx2»Nl 

R1»0,5*(R( J,K)*R( J*1,K)  ) 

R2«U,5*(EX(j,K)*EX(J*l.Kn 

FlaO,5*(r{ J,K)*F(J*1,K)> 

F2»0,5*(EF<j,K)*EF(J*l,K)  ) 

R3iiOi5«(R(  J,K)*R(JiK*l)  ) 

R4bU,5*(EX(J,K)*EX(J,K*1)) 

T380i5*(T( J«K)*T( JiK*l)) 

T4bO,5*(ET(J,K)*ETU.K*1)) 

CFlgO,5*(CF(J)*CF(J*i)) 

SPsO,5*HS 

SVsSP-1.0 

X1»SP*(EX(JjK)-R(J,K))*R(J,K) 

X2«HA*X1*CF(J) 

F6Br(J,K)*SV-SP*EFU,K) 
T2  =  T(J,K)*SV.SP*ETU.K) 

Al=(HC*Xl*X2*(l,0*(F6*F6)/(Xl*Xl)*«t2*T2)/<X?*X2))W(E>{J,K)-RCj,K 
X)\ 
B2(l)«Bl(l) 

Bl<i,)s(RA*Rl*CFl)*(R2-Rl)/Rl 
C2{|,,J)sCl<l.J) 

CKli  J)s(R3*(R4-R3n/(RA*R3*CF(J)  ) 
D2(t)sDl(l) 
XlsSP*(R2-Rl)*Rl 

Dl(l,)a(HD*(RA*Xl*CFl)*(ri*SV-SP*F2)  )/Xl 
E2(],,J)8E1(1,J) 

Xli:SP*(R4«R3)*R3 
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El(;iJ)=(HD*Xl*(SV*T3*SP*T4))/(RA*Xi*cr(J)) 

GAM?|Al*o,5*(Pl(l)*R2(l)*ri(l.J)*C2(l.  J))*Dl(i)-D?(t)*Ei  (li  J)-E2(l, 
IJ) 

WPs(Al*G(2,«,K)*0,5*(Bl<l)*G(l.J*lJK)*B2{l)*S(liJ.i,K)*Cl(l.J)* 
lG(ljJ,K4.1)*C?(l,J)*G(l,J,K-l))*Dl(l)*G(2.J*l'.K)-P2ll)*f{2#J'l#»<)* 
2EX(liJ)*G(2,J»K*l)-E2(l,J)*G(2,j'.  K«1))/GAM 

UP=BEO*G(liw.K)*RELI*WP 

UAbUP»G(1, J,K) 

UDsABStUA) 

UEsAMAXKUEjUD) 

G(li J,K)3UP 

WPs|Al*A(2i«,K)*0,5*(Bl(l)*A{l.J*i;K)fB2(l)»*(l,J»J:iK)*Cl(l#J)* 
lA{ljJ,K*l)*C?(l,J)*A(l.J,K-l))*Dl(l)*A(2iJ*l'.K)»C2tl)**(2»J'l.K)* 
2El<liJ)*A(2jJ.K*l)-E2(liJ)*A(2,J,K"l))/GAM 

UP3HE0*A(l,w.K)*RELI*WP 

UAsgPoAd*  J»K) 

UDsABS(UA) 

UE=AMAXl(UEsUD) 

A(l| J|K)sUP 

DO    14     I«2»NI1 

Si(  l-l)*HS 

SP«S*0.5*HS 

SSsS-1,0 

ST«SPsl,0 

A23A1 

Xl5bP*(EX(J,K)-R(J,K))*RtJiK) 

X2aHA*)(l*CF(  J) 

F'6sr(  J,K)*ST.SP*EFU,K) 

T2sf (J,K)*ST»SP*ET(J,K) 

A1»(HC*X1*X2*(1,0*<F6*F6)/(XI*X1)*<T2*T2)/(X?*X2))W(E>(J|K)- 
IP{J|K)) 

B2(  \)»BH  I) 

XlaiS*(R2-Rl)*Rl 

Bl(|)«(RA*Xl*CFl       )*(R2.R1)/X1 

C2(|i J)8C1( I.J) 

XlsS*(R4-R3)*R3 

Cl(|iJ)sXl*(R4-R3)/(RA*Xl*CF(J)) 

D2{|)»D1(I) 

X1"SP*(R2-R1>*R1 

Dl(  I  )3(HD*<RA*X1*CF1)*(F1*ST-SP*F2) )/Xl 

E2{I,J)sEl(I.J) 

Xl=SP*(R4«R3)*R3 

El(  |iJ)s(HD*Xl*(T3*ST-SP*T4)  >/(RA*>'l*CF(J)  ) 

GAM§A1*a2*B1(  I  )*B2(  n*Cl(  !•  J)*C2M«  J)*D1(  I  )»P2(  I  JbBK  I-1)*D2(  I"!)* 
lEK  l#J)»E2(  I,  J)-Ein-1,J)*E2(  I=1',J) 

WPs(A1*G(I*1.J.K)*A2*G(I-1iJ»K)*b1{I)*G(I,J*1,K)*B3(I)*G(IiJ-1iK)* 
ICKI,  J)*G(  I,  J,K>1)*C2(I,.J)*G{  I,J',K'1)*D1(I)*G(  I*i,J*l,i«)5Dl(  Nl>* 
2G(  I.l,  J*liK)-D2(  I)*Gn*l.J»l»K)4D2<  I-1)*G(  I-i,  j-1,  K)*E1  (  Ii  J)* 
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3G(I#l,J,K*l)-E?(I,J)*GM*liJiK-l) 
4*Q( 1-1, J,K-1)  )/GAM 

UPbHEO*G( liw.K)*RELI*WP 

UAsUP-G(  li JiK) 

U0»A8S(UA) 

UE=4MAX1(UE,UD) 

G(l| J,K)8UP 

WPs1A1*a(I*1,J,K)*A2*A(I-1jJiK)*B1{1) 
ICK  Ii  J)*A(  I,  J,K*1)*C2{  1#  J)*A{  I,  j'.K-l) 
2An"l#  J*1»K)-D2(  I)*A(I*1.J-1iK)*d2(  I 
3A(I*l,J,K*l>-E2(l,J)*A(I*liJiK=l)-El( 
4*A( l-l, J,K-1)  )/GAM 

UPa8E0*A< Iju.K)*RELI*WP 

UA=UP-A( 1, J,K) 

UD«ABS(UA)  -^ 

UE»AMAX1{UE,UD) 

A(l,J,K)aUP 
14   CONTINUE 

BiJ{Ml)8Bl(NI  ) 

Bi(Ml )B(RA*R2*CF1)*(R2-R1)/R2 

C2(M,  J)xCl(NJl,  J) 

CKNii  J)sR4*(R4-R3)/(RA*R4*Cr(  J)  ) 

GAMgAl*0.5*(Bl(Nl )*B2(NI)*C1(NI, J)*C2 
1E2(NI1, J)-El(Ntl, J) 

WP=|A1*G(NI1, J,K)*0.5*(Bl(NI)*G(Nli J* 
ICKM,  J)*G(M,  J,K*1)*C2(NI.  J)*G(Nr»  J. 
2Dl(Nll>*G(NIl,J*l,K)*E2(Nll#J)*G(Nll, 
3)/GAM 

UPbRE0*G(N1, J,K)*RFLI*WP 
UA«UP-G(NI,J,K) 
UDsABSCUA) 
UEsAMAVKUE.UD) 
G{N|,J,K)sUP 

WPa(Al*A(Nll, J,K>*0,5*(Bl(NI)*A(Nlij* 
ICKNI,  J)*A<M,  J,K*l)*C2(Nli  J)*A(NI»  J. 

2D1(NI1>*A(NI1, J*l,K)*E2(Nll,J)*A(Nll, 
3)/GAM 

UP««EO*A<N!, J,K)*RELI*WP 

UUsABS(UA) 

UEsAMAXKUE.UD) 

A(N|, J,K)«UF 
13   CONTINUE 
12   CONTINUE 

PERIODICITY  CONDITIONS  FOR  POTENTIALS 

DO  17  Isl,Nl 

DO  18  K*liN2 

A(I|1,K)«A( I.Nl.K) 

A(I|N2,K)»A( Ii2.K) 


El(I-l,J)*e(I-i,J,K*i)*E2(I-l,J> 


♦A<I,J*1,K)*B2(l)*A(I.J-liK)* 
♦  D1(I)*A(I*1,J*1.K)  =  D1(U1)* 
.1)*A(  Ni,  J.1,K)*E1  (Ii  J)* 
I-liJ)*A(I-l, J,K*l )*E2(I.1,J) 


(NI,J))#D2(Nl!)-Dl (Nil)* 

l.K)*B2(NI>*GlNI,il-l,K)* 

K-l>)*D?(MntG(Nfl,J»l,K)- 

J.K-l)-Fl(N!l^J)*f(NIl,J,K*i) 


1,K)*B2<NI)*aJNI.iL-1i"<)* 
K<»l))*D2{KIt)*A(NfliJ  =  liK)- 
JiK-1)«E1<N!i7j)**(NI1,J,K*i) 


IN  VACOUK  REGION 
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19 
17 


320 


G(ljl,K)sG(I.Nl,K)-1.0 

G(  I|N2,K)=G{ Ii2,K)*1.0 
18   CONTINUE 

DO  19  Jsl,N2 

A(  I|  J,1)*A(  I,  J.ND-l.O 

A(  Ii J,N2)*A(  I» J,2)*1.0 

G(  Ij Jil)sG{ I, J.Nl) 

G(  I, J,N2)sG( 1, J,2) 

CONTINUE 

CONTINUE 

NTE«=NTER*1 

ITEH»ITER*1 

CHECK  ERROR  IN  THE  SOLUTION 

IFCUE.GE.ERRI )G0  TO  200 

UJNgO.O 

SV"0,5*HS 

DO  400  J  =  2,M 

Rl»0,5*(R( J,2)*R(JiNl)) 

R2sQ,5*(T( J,2)*T( J.Nl)  ) 

EKli  J)  =  -(SV*SV*Rl*R2*HD)/CRA*SV*Rl*CF(  J)  ) 
400   CONTINUE 

ITEHATE  SOLUTION  AT  SINQULAR  POINTS 

DO    401    Ks2,M 

SUMX=0.0 

SUM^'O.O 

Rl=0i5*(R(2iK)*R(Nl.K)) 

R2«0i5»(F(2,K)*F(Ni,K)  ) 

R3sO,5«(CF(2)*cr(Nl)  ) 

Dl(X)a.(HD*R2*(RA*SV*Ri*R3) )/Rl 

DO    402    J>2,M 

Xl5BA*SV*R(w.K)*CF( J) 

Al«MC*SV*Xl*(l,0*(F(J,K)*F(JiK))/(P(J.K)*R(J,Kn*(SV*SV-*T(J,K)* 
IT(J,K))/(X1*X1)  ) 

D2(l)aDl(l) 

Rl«0i5*(R( J»K)*R( J*1,K)  > 

Dl(l)a-(HD*0,5*(F(J,K)*r(J*l»K))*(RA*SV*Rl*0T5*(CF|j)*rF(J*l))))/ 
iRi 

E2(liJ)sEl<l, J) 

R4  =  0,5*(R( JiK)*R( J,K*1)  ) 

El{liJ)s-(HC*SV*SV*0.5*(T(J|K)*T(J'K*l))*R3>;(RA*SW*R3*CF(J)) 

SUMlsSUMl*Al*U(2.  JiK)*Dl(l)*U(2,  J*l,K)=D2(l)«U(2i  Jil.KWEKl.  J)* 
lU(2iJ,K*l)-E?(l,J)*U(2,J,K-l) 

SUM2«SUM2*A1*D1(1)»D2(1)*E1<1, J>-E2(1, J) 
40^      CONTINUE 

AASSUM1/SUM2 

UAbAA«u(1#2,K) 

UU=ABS(UA) 

UINIAMAXKUC.UIN) 
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DO  403  J«2,M 
404   U(1|J,K)»AA 
401   CONflNUE 

ITEWATE  SOLUTION  FOR  POTENTIAL  U  I^  PLASMA  REGION 

DO  402  I«1.M 

SV»ll-0,5>*HS 

DO  301  J«2,M 

Rl8U,5*(R( Js?)*R( J.Nl)  ) 

R280,5*{T(Ji?)*T(J.Nl)) 

Sa( l-l)*HS 

CX( li J)s(S*Rl*Rl)/(RA*S*Rl*Cr( J)  ) 

EX(i#J)s-<SV*SV*Rl*R2*HD)/(RA*SV*Rl*cr(J>) 

30X   CONTINUE 
302   CONTINUE 

DO  412  Kb2,M 
R3«»0i5*(Cr(2)*CF(Nl)  ) 
Rl'«0.5*(R(2»K)*R(NliK)  ) 
R2«»0i5*(F(2iK)*F(Nl,K)) 

R4lRl*R3 
R5«HD*R2 
DO  !»20  I«1,M 
Sf ( J-1)*HS 
Bl(  l)«(RA*S*R4)/S 

520   CONTINUE 

DO  S21  I«1,M 

Sf ( l-l)*HS 

SPsS*0.5*HS 

Dl< I >«-(R5*(RA*SP*R4))/Rl 

521  CONTINUE 

DO    413    Jb2,M 

Rl*0i5*(RUiK)*R(J*l,K)) 

Fl«0,5*(F(JjK)*F(J*l,K)) 

R4aUi5*(R(JiK)*R(JiK*l)) 

T330,5*(T(JbK>*T(J.K*1)) 

CFH0,5*(CF(  J)*CF(J*1)) 

SP8Q,5*HS 

F28P(J,K)*F{ J,K) 

T2bT(J,K)*T(J,K) 

RR2iR( J,K)*R( J.K) 

XlsHA*SP*R(J,K)*CF(J) 

AlaHC*SP*Xl*(l,0*F2/RR2*(SP*Sp*T2)^(Xl*Xl)) 

D2(1)«D1(1) 
Dl(l)i-(HD*Fl*<RA*SP*Rl*erl) )/Rl 

E2(4,J)«E1(1.J) 

EKli  J)s-(HL*SP*SP*T3*R3)/<RA*SP*R3*CF(J)  ) 

DO    414    I»2,M1 

Si(l-1)*HS 

SP8S*0,5*HS 
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SS»SP*SP 

Xi3HA*SP*R(»<,K)*CF(  J) 

Al8HC*SP*Xl*(1.0*r2/RR2*(SS*T2)/(Xl*Xl)) 

B2(|)nBl<  I) 

Bl(  |)»(RA*S*Rl*Cri)/S 

C2(  |,J)sCl(  I.J) 

Ci(  |i J)s(S*R3*R3)/(RA*S*R3*Cr( J) ) 

n2<  PiDKl) 

DK  I  )s.(HD*Fl*(RA*SP*Rl*CFl) )/Rl 

E2(liJ)«El(l,J) 

EK  Ii  J)»-(HC*SS*T3*R3)/(RA*SP*R3*C''(  J>) 

GAMiAl*A2*Bl(  I  )*B2(  n*Cl(  I»J>*C2(  I»  J)*D1<  I)»C2(  D.BK  I.l)*n2(  I"l)* 
1E1(|,J)-E2(I,J),E1(I-1,J)*E2{I=1.J) 

WPsilAl*u(I*ltJ.K)*A2*UM-l»JiK)*Bl(l)*U(IiJ*i,K)*B2(I)*UniJ-l»K)* 
lCl(J,J)*U(I»J.K*l)*C2{I,J)*U{r.J.K-l)*Dl(I)*0(l*i,J*l.l«)=Dl(I»l)* 

2U{I«l,J*l,K)-D?(l)*U(I*l.J«l#K)*D2(I-l)*U{l-l,j-l,K)*El(I.J)* 

3U(  1*1, J,K*1)-E2(  I, J)*U(I*l»JiKcl)»fel(I-li J)*C(I-li J,K*l )*E2( I-l, J) 

4*U( (-1, J.K-1) )/GAM 

UPeHAO*U(  liW.K)*RULA*WP 

UAsgP»U( 1/ J»K) 

UU3ABS(UA) 

UlNBAMAXKUCUIN) 

UniJiK)»UP 
314      CONTINUE 

B2(NI )>B1(NI  ) 

Bi(Nl)«RA*Rl*Cri 

C2(NI, J)aCl(Nl, J) 

CKM,  J)  =  (R3*R3)/(RA*R3*CF(J)  ) 

GAMiAl*o,5*(Rl(NI)*B2(Nl)*Cl{Nl, J)*C?(NI,J))«D2(Nli)»Di (Nil)* 
XE2(NI1,J)-E1(NI1, J) 

WP«|Al*u(Nll,  J,K)*0.5*(Bl<NI)*u(Nli  J*1,K)*B2(NI)«UINML-1«'<)* 
ICKNli  J>*U(Mi  J.K*l)*C2(Nli  J)*U(NI»  J,K  =  1)  )  ♦D?  {  M  1)  *U  <  N  1 1 .  J  =  l ,  K  ) - 
2Dl(Nll)*U(NIl,J*l,K)*E2(NlliJ)*U(Nll,J,K-l)*'El(N!l7j)«l(NIl.J|K4.i) 
3)/GAM 

UPsHAO*U(NI, J|K)*RULA*WP 

UAbUP.U(NI,„,K) 

UO»ABS(UA) 

UlNsAMAXKUlNiUD) 

U(N|,J,K)«UF 
313   CONTINUE 
312   CONTINUE 

PERIODICITY  CONDITIONS  FOR  POTENTIAL  IN  THE  PLASMA  REGTON 

DO  317  Isl.M 

DO  318  Ksl.N? 

U{I,1,K)»U{I,N1,K) 

U(I|N2,K)«U(I».2,K) 
318   CONTINUE 
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DO  J19  J«1,N? 

U(Ii  J,1)»U(I,  J.ND^l.O 

U( I| J,N2)=U( 1; J,2)*1.0 
319   CONTINUE 
31?   CONTINUE 

iNEHslNER*! 
C         CHECK  ERROR  !N  THE  SOLUTION 

IF(UIN.GE,UE)G0  TO  320 

BMAXaO.O 

SUMiO.O 
C         COMPUTE  VOLUME  OF    PLASMA  REGION 

DO  SI  Ja2iNl 

F"l3iCF(  J)*TER 

DO  SO  K82#N1 

AUXiR( J,K)*R( J,K)*(RA2«R<JiK)*FI3) 
SUMiSUM*AUX 

50  continue 
5x  continue 
volbh*h*sum 
c      compute  pressure  in  the  plasma 

gra8(2,0*con)/(vol**gas) 
c      compute  currents  01,02.03 

CALL  rLAX(rLUl,FLU2,FLU3.Ql.Q2.Q3) 
C         COMPUTE  ERROR  IN  THE  BOUNDARY  AND  l-SE  STEEPEST  DPS8ENT  METHOD  TO 
C         UPDATE  BOUNDARY  POSITION 

CALL  GRAD(BKAX.PR01,Q1.Q?,Q3) 

PR0l'PR01/<NJ*NJ) 

BET»GRA/PR01 

PRINT  21,  ITER,V0L,RRA,UE.UlN,BMAX,f^l,Q2,03.BET.ENER 

21  FORMAT  (^  ^,2X,  M,6X,F8;3.4X,F6;3'.3X,E10,4,6X;e10,4;6X,F10,4,3Xi 
1F5,1»3X,F5,1,3X,F5,1.3X,F5,3,3X,E11.6) 

IF(NTER.GE,100)GO  to  23 
C         COMPUTE  ERROR  TO  BE  ALLOWED  IN  SOL^'TION  IN  T|iE  INTERIOP,AS  A 
C         FUNCTION  OF  THF  ERROR  IN  THE  BOUNDARY 

BMN8FAC*BMAX/Ql 

ERRJaAMINK  ARRI.BMN) 

CALL  SECOND(TIM) 
C         CHECK  TIME  ELLAPSED.AND  IF  GREATER  THAN  TOM, TERMINATE  THE 
C         COMPUTATION 

IF(TIM,GE,TCM)  GO  TO  1570 
C         CHECK  BOUNDARY  ERROR 

IF(BMAX,GE,EPRR)RO  TO  IQO 

PRINT  22,  (  (R(  J,K),  Js2,Nl,N3),K  =  2'.  Nl,N3) 

22  FORMAT(^  ^,13F7,3) 
IFdND.LE.O)  STOP 
IF{BEF,LE,0,0)STOP 
GO  TO  25 

23  CONTINUE 
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PRINT  32 
32  fORHI^J{^l/,/FR^E    SURFACE?*) 

PRINT  22i  <  {R(  JiK),Kr2,Nl,N3)|  Js2'.  M) 

PRINT    35,AH,A12.A22.A33 
35      FORHMt/    ?!,3X,4E15,4) 

NTE8»0 

GO  TO  100 
2;>   CONTINUE 
C         REFJNE  THE  KFSH.AND  COMPUTE  INITIAL  SOLUTION  BY  FXf RAPf LAT !  ON 
C         OF  SOLUTION  FOR  THE  PREVIOUS  MESH 

HlO,5*H 

HS«g,5*HS 

HAaO,5/HS 

HB=0i5/H 

HCs(H*H)/(HS*HS) 

HDbM/(2,0*HS) 

DO    800     II»1«NI 

I*N|-I1*1 

DO    800    J52,M 

DO    800    K«2,M 

A(2«I-l,J#K)eA( I.J.K) 

U(2tl-1,  J#K)t:U(  I.J.K) 
800      G(2«N1,  JiK)bG(I.  J,K) 

DO    801    1*1, Ml 

DO    801    js2,M 

DO    801    K«2,M 

A(2«I,J,K)»0,5*(A(2*1-1,J,K)*A(2*I*1;J.K)) 

G(2«l, J,K>sO,5*(G( 2*1.1, J,K)*G (2*1*1. JiK)) 
801    U(2*Ii J,K)s0.5*(U(2*I.l, J,K)*U(2*I*1.J.K)) 

NJb2*ni=1 

N|HNI,1 

DO    802    JJ»2,Nl 

JlNl*2«JJ 

DO    802    K«2,M 

R{2«J-2,K)eR( J,K) 

DO    802    1=1, M 

A(Ij2*J-2iK)eA( I, J,K) 

G(Ij2*J-2#K)aG( I, J.K) 
fl02    U( I,2*J-2iK)8U( I.J.K) 

DO    803    J«2,NJ 

DO    803    K«2,M 

R(2tJ-l,K>»0.5*(R{2*J,K)*R(2*Js2.K'  ) 

DO    803    I«1,M 

A(Ii2*J-liK)s0.5*<A(I,2*J|K)*A(I^2*J-2iK)) 
G(I|2*J»liK)sO,5*(r,(I,2*J#K)*G(I.2*J"2iK)) 
80J      U(I|2*J  =  l,K)r0,5*(U<  i,2*J,K)*U(  r,2*J-2iK)  ) 
DO    810    K«2,M  , 
R(2«N1.i,K)50,5*(R(2.K)*R(2*N1s2,K5  ) 
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DO    810    I«1,M 

A(I,2*Nl-l,K)30,5*<Aei.2.K)*A(I,2*^l«2.K)) 

G{  I|2*Ni-l,K)sO,5*(G(  r,2.K)*G( l,2*Ni,2iK)*l,0) 

810  Uni2*Nl-liK)so,5*(un.2.K)*U(l,2*Nl-2.K)) 
NN«2*NJ 

NNHNN*1 

DO    «04    KKB2|Nl 

K»Ni*2-KK 

DO    804    J»2,NN1 

R( J|2*Kc2)=R( J.K) 

DO    804    Isl,M 

A(  li J,2*K-2)xA( I, J,K) 

G(  Ii Ji2*K-2)=G{ 1  I JiK) 
804      U{I,  J,2*K-2)*UniJ.K) 

DO    805    Ks2,Nj 

DO    805    J  =  2,NN1  "^ 

R( Ji2*K-l)«0,5*(R( J,2*K)*R(Ji2*K-2) ) 

DO    805    I«1,M 

A(IiJ,2*K-l)»0,5*(A(I,J,?*K)*A(I,Ji2*Ko2)) 

G(I|Jj2*K-l)  =  0,5*(G(  I, J,2*K)*G(  I. J»2*K-2) ) 
80t>      U(I,J,2*K-l)  =  Q,5*(U(l,J,2*K)*U(r,J*2*K  =  2)) 

DO    811    Js2,NMl 

R( J|2*N1-1);0,5*(R( J,2)*R( J#2*N1=2) ) 

DO    811    Isl,M 

A(I|J,2*Nl-l)sO.5*(A(I,J.2)*A(I,j,2*Nl=2)*l,0) 

G(IiJ,2*Nl«l)=0,5*(G(r.J,2)*G(l,j,2*Nl=2)) 

811  U(  I| J,2*N1-1)«0.5*(U(  r,Ji2)*U{ I, j,2*Nlc2)*l,0) 
NJa2*Nj 


N1bNJ*1 

N2«NJ*2 

CALL  SARF 

DO  807  Isl.M 

DO  808  K»1,N? 

A(  lil,K)«A( I.Nl.K) 

U(1|1,K)«U(I.N1,K) 

Gnil|K)»G(I,Nl,K).l 

.0 

A(I|N2,K)«A(I,2,K> 

UniN2,K)"U(  I.2,K) 

808 

G{l|N2,K)"G(I,2,K)*l 
DO  809  J»1#N2 

.0 

A(I|J,l)sA(I.J,Nl>.l 

.0 

U(I|J,1)bU(  I.J.ND-l 

.0 

G(1|J,1)bG< I.J.Nl) 

A(  I| J,N2>»A( I, J,2)*l 

.0 

U(I,J,N2)«U{li J,2)*l 

.0 

809 

G(I|J|N2)»G(I, J,2) 

807 

CONTINUE 
REFi-1.0 
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1570 


N384 

ERRa=ORRB 
ERR|»ORRI 
FACiFEC 

REL|=RELl*0,5*(2.0-RELn 
RULA»RULA*0,5*(2.0«RULA> 

REOsl.O-RELl 
RAOhI.O-RULA 
GO  TO  2000 

CONTINUE 
REWIND  2 
WRITE  RESULTS 


W«ITE{2> 

WRITE(2) 

WRITE<?) 

WRITE<2) 

WRITE(2) 

WR1TE<2) 

STOP 

END 


((  (A 


IN  TAPEl 
( Ii J.K).Ji 


il.N2),K.l'.N2).Ial,NI) 
(({U(  1,  J#K).J«l.N2),Kal'.NH),  I=1,NI) 
(((R(  I,  J,K),Jo1.N2),Ks1',n2)'.  IsliNl) 
(  (R( JiK), J81,N2)|K  =  1#N2) 
( (EX(J,K), Jsl,N2)iKal,N2) 
FLLl#rLU2,FLU3,C0N 


C 
C 
C 


SUBBOUTINE  ASYM( FLU3 . FLU2i FLUl. ASY^ ) 

COMPUTES  THE  VALUES  OF  THE  FlXED  FLyXES  FLUl^FLU?, PLU3 . IN  TERMS 

OF  THE  VALUES  OF  THE  INITIAL  CURRENTS .  COMPUTES  A^ULLY  SYMMETRIC 

SOLUTION  WHICH  WILL  BE  USED  AS  INITIAL  GUESS  FOR  GENERAL  SOLUTION 

COMMON  GG(9,52,52),GX(9,52#52),GY(9,52.52),Gi(52.52).Q?{52i52)# 
1G3(52,52)»G(9,52),A(9,52).B(9,52).P(9,52),CF(  52  )  ,  P*  <  52  »  ,  R  (  52  ) , 
2EX(5»2),EF(52)iRA.RO,Rr.NIiNJ,PETE,^AFI,BETE,GRA,N,^S.VfL.ENER, 
3RELb 

DlMbNSION  CP(52) ,CV(52) 

PEF|Bl,0 
305   FORMATC^  ^,3x,3Fl5,3) 

EHR|«0, 00001 

RELI=1.9 

ICOsl 

ACOil.O 

ERRa=ASYE 

ARR|=ERRI 

BMAXsl.O 

Pl=4, 1415926535898 

F|TSPETE/(2,0*PI  ) 

FTP|8ETE/(2,0*PI) 

FTViPETE/(2.0*PI) 

PPsl,0/{2,0*Pl) 

H»(2,0*PI  )/<NJ-0.0) 
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HS»1,0/(NI«1.0) 

NlsNJ*l 

N2sNJ*2 

HAsO,5/HS 

HB"0i5/H 

HCa(H*H)/(HS*HS) 

HU=H/(2.0*HS) 

RtO«l,0-RELl 

NIllNI-l 

VOLI|2,0*PI*PI*RI*RI*RA 

GASi5, 0/3,0 

CONiO,5*GRA*(VOL**GAS) 

RA2JiO,5*RA 

CALL  ELI 
DO  400  js2,M 
300   EF(d)aHB*(EX( J*1)-FX(J-1)  ) 
DO  i    J*1,N2 
DO  I  lal.NI 

1  G{ 1| J)«PP*PEFI*( J-2)*H 
DO  2  J»liN2 

2  CF(J)»C0S((J-2)*H) 

SUMI»0.0 

suMaso.o 

DO    501    Js2,Ml 
SEMl'O.O 
SEM2=0.0 
DO    502    Isl.Ml 
S»(  |»1)*HS*0,5*HS 
Rl5S*R(J) 

R2s(EX( J)-R( J) )*S*R( J) 
SeMl.»SEMl*Rl/(RA*Rl*CF(J)) 
SEM2»SEM2*R2/(RA*R?*CF(J)) 

502   CONTINUE 

SUMl,aSEMl*R(  J)*SUM1 
SUM2=SEM2*(EX( J)-R(J) )*SUM2 

50i   CONflNUE 

SUMiaSUMl*H*HS 

SUM2sSUM2*H*HS 

FL,U2«FTV*SUK? 

Fl,UlsFTP*SUf'l 

ITEHsQ 

NTEH»0 
100    R(l|sR(Nl) 

R(N2)»R(2) 

EF(l)aEF(Nl) 

EF(N2)«EF(2) 

DO    4    Js2.Nl 
3         F( J|sHB*(R( J*1)«R( J"l)  ) 
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5 

4 


6 

200 


F{ 

DO 

Ri 
Fi 

R2 

r2 

CF 
CV 
CP 
DO 

S3 

SP 
XI 
X2 
A( 
XI 
B« 
XI 
X2 
D( 
CO 
CO 
CV 
CV 
CP 
CP 
DO 
B( 
B( 
A( 
A( 
D( 
D( 
UE 
DO 
GA 
WP 

1G( 
WP 
UP 
UA 
UU 
UE 
G( 
DO 
GA 

ID( 


H=F( 

N2)»r 

4  J« 
=»0i5* 
=  0.5* 
sO,5* 
=  0.5* 
1«0,5 
<J)'( 

&  I« 
(I'D 


3bs*o,; 

L=SP*<E 

'IiJ)a' 

=  S*(R 

I|J)« 

sSP*{ 

«P1*( 

ljJ)  = 

IN?  INUI 

IN?  INU 

(l)aC 

(N2)« 

(l)»Cl 

(N2)s 

6    Is 

I|N2)! 

I|1)B 

I|N2) 
III)* 
I|N2) 
=  0.0 
i  ?  J« 
iM8A(l 

■  Adi 

2|J*1 

aWP/G 

"WEO*, 

k»UP-Gl 

)  =  ABS(l 

^bAMAX" 

llJ)  = 

I    8    Is 

iMrA(  1 

I»li  J 


Nl) 

(2) 

2.N1 

(R(J) 

(F(J) 

(EX(w 

(EF(j 

*(CF( 

R(J>* 

R(  J)* 

l.NI 

*MS 

5*HS 

EX(  J) 

*(SP« 

(HC*X 

2-Rl) 
(RA*X 

R2-R1 
SP-1, 
(HD*X 

E 

E 

V(N1) 

CV{2) 

P(N1) 

CP{2) 

l.NI 

B(I#N 
sB(  If 
AdiN 
sA(  I( 
D(IiN 
sD(  I, 


♦R( J*l) ) 

♦  r( J*l)  ) 

)*EX( J*l)  ) 

)*FF( J*l)) 

J)*CF( J*l)  ) 

(EX(J)-R(J)>)/(RA*R(J)*CF(J)) 

P( J) )/(RA*R( J)*CF( J)  ) 


-R( J) )*R( J) 

1.0)-SP*EF(J) 

l*(RA*Xl*CF(J))*(1.0*(X2*X2)/(Xl*Xl)ni(EXtj)*R(  J)) 

♦  Rl 

1*CF1)*(R2-R1 )/Xl 

)*R1 

0)-SP*F2 

?*(RA*X1*CF1) )/Xl 


1) 
?> 
1) 
?) 
J  ) 
?) 


2,N1 

.J)*0 

J)*G( 

)-D{l 

AM 

G(li J)*RELI*WP 

(1>  J) 

UA) 

1(UE,UD) 

UP 

2»NI1 

,J)*A(I-l,J)*B{l,J)*Bn.J-l'*D(IiJ)-Dn-liJ»-D(!,J-l)* 

1) 


.5*(B(1,J)*B(1.J-1))*D(1.J)»D(1,J«1) 

2iJ)*0,5*(B(liJ)*G(l'.  J*l)*B(l»J*l)*GflTj  =  in*D(l.j)* 
,  J-1)*G(2, J-1) 
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61 


62 

60 


WP3A(I,J)*G(Ifl,J)*A(I-l,J)*G(l.i,J)*B(IiJ)*S(l,J*i)*BniJ»l)* 
iGdiJ-D^DCUJJ^Gd+l.J^D-DCNl'.  J)*G(I-l,J*n^D<l7J  =  l»*G(!*l,J-l) 
2*D( 1-1, J-1)*G( I«l, J-1) 

WPsWP/GAM 

UPsREO*G( !• J)*RELI*WP 

UAsgPcG< li J) 

UPaABS(UA) 

UE*AMAX1(UE,()D) 

G{I|J)*UP 

CONf  INUE 

GAM?A(NI-1,  j)  +  0.5*(8(Nl,  J)+B(NI,  J-1)  )*D(NI-1',Jt1),B(NI-1,J) 
WPsA(NI-l#J)*G(Nl-l,J)*0.5*(B(Nr.j)*G(NI,J*l)+B(Ni;j=l»*G<NI,J-l)) 
l*D(NUl,  J-l)*G(NIli  J-l)-n(NI-l,  J)*G{NI1,  J*l) 
WPsWP/GAM 

UP5SE0*G(NI, J)*REL1*WP 
UA«UP»G(NI,J) 
U13sABS(UA> 
UE  =  AMAXl(UE,llD) 
G(N|,J)sUP 
CONf INUE 

DO  9  lal.NI  ' , 

G(Ij1)bG(  I,M)-PEFI 
Q(I|N2)=G(1,?)*PEFI 

1TER»1TER*1 

If (UE.GE.ERRI )G0  TO  200 

NTE«aNTER*l 

ENERsO.O 

DO  *0  Js2,Nl 

DO  61  IsliNIl 

X1»U( 1*1, J)in( I, J) 

X2aG( 1*1, J*l)-G( I, J) 

X3«li(I*l,J)»G(I,  J  +  1) 

ENEH»ENER*A(  I  ,  J  )  *X1  *X1*D(  I  •  J)  *  (  XS^^'S-XS^XS  ) 

XlsQd,  J*1>»G(1,  J) 
ENEH=ENER*B(li J)*X1*X1*0.5 

X1=U(NI, J*l)-G{NIi J) 
ENEWsENER*B(tvJl,  J)*xi*Xl*n,5 

no  62  Is2#Ml 

Xlslid.J^D-GC  I,  J) 

ENEHsENER*B( I , J)*X1*X1 

CONTINUE 

ENESs(ENER*HS)/H 

ANE«eENER*2in*PI 

BMAX'O.O 

5UM|0,0 

IF(ACO,LE,0,n)  GO  TO  500 

ACO»»1,0 

F|,Ui  =  ANER*PAFl 
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500       Al,A»rLU3/ANiER 

ENEKsO,5*ANER*ALA*ALA*0,75*GRa*VOL 
CALL    FLUX{FLUl,FLU2,rTV,FTP) 
ENEfil»ENER*PI*(FTP*FLUl*rTV*FUU2> 
DO    10    Js2iNl 
AUXiR(J)*R(J)*{RA2*(R(J)*CF(J))/3.0) 

10  SOM(|SUM*AUX 
VOLi2,0*Pl*h-*SUM 
GRAS(2.0*CON)/(VOL**GAS> 
DO    11    Js2,Nl 

A2  =  1.,0/(RA*R(  J)*CF(  J)  ) 

XlaG(2, J-1)«R(1, J-l) 

X2=G(2, J)-G(l» J) 

X3ali(2,  J)-G(2»  J-1) 

X4sG{l, J)-G(l , J-1) 

X5«G(2, J)-G(l, J-1) 

X6«lj!(l,  J)-G(?,  J.l) 

X7««G(2,  J)"G(li  J) 

X8sU(2, J*1)*G<1, J*l) 

X9=G(2, J*l)-n(2. J) 

XlOgGd,  J*l)-G(l,  J) 

XllfiG(2, J*1)-G(1,J) 

Xl2nG(l, J*l)-G(2, J) 

SAMgO,5*(0,5«-(A(l,  J-l)#Xl*Xl*A(l'.j'*x2*X2)*0l5*(B(f,  J»1  )*X4*X4* 
lB(ljJ=l)*X4*x4)*D(l,J.l)*(X5*X5-X6*X6)*0,5*(A(l,J)*X7*y7*A<l.J*l)* 
2X8*X8)4.  B<1.  J)*XlO*xiO    *D(l'.  J)*{Xtl*X11  =  Xl2*X12)  ) 

R1=H*H*R( J)*(RA*R{J)*CF( J) ) 

SAMiSAM/(Rl*(EX( J).R(J) ) ) 

VAMi|H*H*(2,0*CV{  J)*CV(  J-D^CVC  j*i)  )*0,25 

PAMiH*H*(2,0*CP( J)*CP( Jpl)*CP( J4i) )*0,25 

VAMiiVAM/(Rl*(EX(  J).R(J)  )  ) 

PAMliPAM/(Rl*R(  J)  ) 

SAM|SAM*ALA*ALA 

VAM?VAM*FTV*rTV 

PAMiPAM*FTP«FTP 

TAMiSAM*VAM 

PRAiPAM*GRA 

D1F>|TAM-PRA 

DAFsDlF/TAM 

CC»ABS(DAF) 

BMAX«AMAX1(BMAX,CC) 

11  R( JJsR(J)-RELB*DIF 

PRINT  12» ITER, VOL, QRA,UE.BHAX,ENER»  ALA,FTV,rTP 

12  FORMAT  (^  ^,3X,  I5i4X,F8.'3.4X,F8.3'.4X,E10,4,4x;e10.4:4X,P12.6,3X, 
13F9,1) 

IF< lTeR.GE,40000)GO  TO  23 
IF{NTER.GE,50)GO  to  27 
BMNS0,1*BMAX/ALA 
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ERRlsAMINK  ARRI.BMN) 

IF(dMAy.GE,ERRR)GO  TO  lOO 

PRINT  22i (R{ J)i J»2,N1,2> 
22   rORMATt^  ^,15r7,3> 

PRIMT  22. (EX(J), Js2,Nl,2) 

PRINT  305»FLlll,rLU?.FLU3 

GO  TO  25 
27    NTEWO 

PRINT  22, (R(J),J=2,N1,2) 

PRINT  22. <EX( J),Js2,Nl,2) 
PRJNT  28 
28   FORMAK^IP!) 

GO  ?0  100 
?.i      CONTINUE 

STOP 
2S   CONTINUE 

PRINT    53 
53         rORMAT(^l;<) 

RETURN 

END 


SUMisO.O 

SUM2=0,0 

NAaNUl 

NlsNJ*l 

DO    I    Js2,Nl 

SEMlaQ.O 

SfcM2»Q.O 
DO    2    1*1, NA 

S»( |"1)*HS*0,5*HS 
RlsS*R{ J) 
P2r— -    -■ 

SEf 

S6( 

CONTINUE 
SUMlaSEMl*R( J)*SUM1 
SUM2»SEM2*(Rn-R( J) )*SUM2 

CONTINUE 

SUMlsSUMl*H*HS 

SUM2sSUM2*H*HS 


■I  |"l)*Hi*U ,5*HS 

lsS*R{ J) 

2s|R0-R( J) )*S*R(J) 

:M3taSEMl*Rl/(RA*Rl*CF(J)) 

:M2sSEM2*R2/(RA*R2*CF(J) ) 

3NTINUE 
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FTV»rLU2/SU^2 

FTPHFLUI/SUM 

RETURN 

END 


C 
C 


SUB«OUTINE  rLAX(rLUl.rLU2iF'LU3.Ql,C2,Q3) 

COMPUTES  VALUES  OF  THE  CURRENTS  Qli0?,Q3#F0R  EACH  POSITION  OF 

THE  FREE  BOlNDARY,lN  TERMS  OF  THE  GjvEN  FLUXES  FLUI , FLL 2 • FLU3 

COMMON  G( 9,52,52), A (9, 52.52 ),U(9;52, 52) •R(52;5?).Fl52,e2).T(52i52) 
1»C1|9,52)#C2(9,52),E1(9,52),E2(9.52).CF(92),P1(52);B2("2),D1(52), 
2D2{S2),RA,RC,RI,NI,NJ,PETEiPEF1,RETe.GRAiH,HS,V0LiENRR.RELB 

COMMON  EX(52,52),EF(92,52).ET(52',52)'.A11iA12;a22,A33 

N«NJ*1 

M«Nl«l 

NlsNJ*l 

HC=|H*H)/(HS*HS) 

HDbM/(2.0*HS) 

AlliO.O 

A12i0.0 

A2280,0 

Aj390,0 

DO  41  Ks2iNl 

DO  i2    Jc2iNl 

Rl=0,5*(R( JiK>*R(J*1,K)> 

R2=Q,5*(EX(w,K)*EX(J*1,K)) 

Flsg,5*(F( J,K)*F(J*1.K)> 

F2sU,5*(EF(j.K)*EF(J*l,K)) 

Ri  =  Q,5*(R( J«K)*R( JiK*l)  ) 

R4aO,5*(EX(^,K)*EX(J.K*l)  ) 

T3sQ,5#{T( J,K)*T( J,K*1)  ) 

T4sU,5*(ET(,j,K)*ET(J.K*l)) 

CFllO,5*(CF( J)*CF(J*1)) 

GSsUd,  j,K) 

A8«A(1,J.K) 

USsUd.J.K) 

G1»Q(2,J,K)»G8 

Al«A(2, J.K)iA6 

Ulsg(2,J,K).u8 

Q2»(j(l,  J*1,K)-G8 

A2sA(l,J*l,K)-A8 

Q3»Q(1, J,K*1)»G8 

A3«A(1,  J.K*1)-.A8 

G48U(2, J*liK)-G8 

A4aA{2, J*liK)rA8 

U4«U(2, J*1.K)-U8 
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Gi>8G( 

U9aU< 
G6aU( 

A6bA( 

A7bA( 
U7»U( 
SP  =  0, 
SVaSP 
XlsSP 
X2sRA 
F6«f  < 
T2  =  T( 
XislH 

XR(J,K 
XlsHA 
Y3eHC 

1T(J,K 
X43|R 
X&ajR 
XlsSP 
XbafH 
Y6>i( 
XXabP 
X7s|H 
Y7«»( 
AlUA 

t(G6»G 
A22SA 

IX7*(A 
A12)|A 

1X7*(A 
A33RA 

DO  J3 

St(I 

SPaS* 

STsSP 
G8aU( 

ASaA( 

U8bU( 

AiBA( 

Ul3U( 
G2bU( 
A28A( 


2,J,K)- 

2,J,K)- 

2, J.K). 

2,J,K4-1 

2,J,K*1 

2,J,K*1 

2,J,K)» 

2, J.K). 

2iJ#K)- 

5*HS 

•1.0 

•(EX(J, 

♦xi*cr{ 

J,K)*SV 
J,K)*SV 
C*X1*X2 

)) 

♦SP*R(w 
*SP*X1* 
) )/(Xl* 

A*Ri*cr 

3*(R4«R 
*(R2-R1 
D*(RA*-X 
HD*ri*( 
*(R4-R3 
D*X1*(S 
HD*SP*S 
11*X3*G 
6«G7*G7 
22*X3*A 
6*A6«A7 
12*X3*A 
6*G6-A7 
33fY3*L 

Is2#M 
1)*HS 
0.5*HS 
1.0 
•1.0 
I, J.K) 
I.JiK) 
I. J.K) 
1*1. J.K 
1*1. J.K 
1*1. J.K 
I.J*1.K 
I.J*1.K 


R(l, J*1,K) 
A(l, J*1,K) 
U(l, J*1,K) 
)-G8 
)-A8 
)-U8 

G(l, J.K+l) 
Ad,  J,K*1) 
U(l, J.K*1) 


K)-R(J,K)  )*R(J.K) 

J) 

-SP*EF(J,K) 

-SP*ET(J,K) 

*(1,0*(F6*F6)/(X1*X1)*<T2*T2)/(X?*X2))W(E>(J.K)- 

,K)*CF(J) 
(1.0*(F(J,K)*F(J,K))/(R(J,K)*R(J;K))*(SP*SF*T(J,K)* 

XI)  ) 

1)*(R2-R1)/R1 

3) )/(RA*R3*CF( J) ) 

)*R1 

1*CF1)*(F1*SV-SP*F2) )/>l 

RA*SP*R1*CF1) )/Rl 

)*R3 

V*T3-SP*T4) >/(RA*Xl*CF<J)) 

P*T3*R3)/(RA*SP*R3*CF( J) ) 

l*Gl*0,5*(X4*G2*G2*X5*G3«G3)*X6*JG4*G4iG5*f5)*X7* 

) 

l*Al*0,5*(X4*A2*A2*X5*A3*A3)*X6*lA4*A4iA5**5)* 

*A7) 

l*Gl*0,5*(X4*A2*O2*X5*A3*G3)*X6*|A4*54,A5*f5)* 

♦  G7) 

1*U1*Y6*(U4*U4<:U5*U5)*^7*(U6*U6-IJ7*U7) 


)-G8 
)-A8 
)-U8 
)«08 
)-A8 
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33 


U2bU 
GJsU 
A33A 
U3sU 
G430 
A4aA 

U43g 

GSsQ 
A;}sA 
Ubsu 
G6sU 
A6sA 
U6  =  U 
G7bU 

A7  =  A 
U7sU 
X1  =  SP 
X28RA 

F65f 
T2  =  T 
X3s|H 

IR(J|K 
XlsWA 
Y3sHC 

1TU|K 
Xlsb* 
X4s(R 
Y4s(R 

XIbS* 
X5  =  X1 
Yi?»|S 

XlcSP 
X6s(H 

Y6»t( 
XlsSP 
X7=|H 
Y7s«( 
AllaA 

IG7*U7 
A229A 

l(A6tA 
A12SA 

l(A6tG 
A33«A 

1(U6«U 
CONf  I 

AlsA( 


«I.J,K*1 
<I.J.K*1 
(I.J.K*1 

{ I*li J*l 

< I*liJ*l 
(Ul.J.K 

(I*li JiK 

<I*liJ.K 

(1*1, J, K 

{1*1. JiK 

( 1*1. J,K 

( 1*1. J, K 

*(EX{J, 

♦X1*CF{ 

(J,K>*ST 

(J,K)*ST 

C*X1*X2 

)) 

*SP*R(u 
*SP*X1* 
))/(Xl* 

(R2-R1) 
A*X1*CF 

A*S*R1* 

(R4-R3) 

*(R4-R3 

*R3*R3) 

*(R2-R1 

D*(RA*X 

HD*F1*( 

♦(R4-R3 

D*X1*(T 

HD*SP*S 

11*X3*G 

) 

22*X3*A 

6-A7*A7 
12*X3*A 
6-A7*G7 
33*Y3*L 
6-U7*U7 
NUE 

NI, J*l, 
NI, J*l, 


)-lJ8 

)-r,8 

)-A8 

)'U8 

,K)-G8 

,K)tA8 

,K),U8 

)-G( I, J*1,K> 

)-A{ I, J*1,K> 

)-U(I,J*l,K) 

*1)«G8 

♦1)-A8 

•fl)-U8 

)-G(  I, J,K*1) 

)-A(  I# JiK*l> 

)-U( I, J,K*1) 

K)-R(J,K) )*R(JiK) 

J ) 

-.SP*EF(  J.K) 

-SP*ET(J,K) 

*(l,0*(F6*F6)/<Xl*Xl)*(T2*T2>/(X2*y2))W(E>(J|K)  = 

,K)*CF( J) 
(1.0*(F(J,K)*F(J.K))/(f^(J,K)*R(j;K))*(SP*SF*T(J,K)* 

Xl)  ) 

*R1 

1)*(R2-R1)/X1 

CF1)/S 

♦  R3 

)/(RA*Xl*CF( J)  ) 

/(RA*S*R3*Cr( J) ) 

)*R1 

1*CF1)*(F1*ST»SP*F2)  >/M 

RA*SP*R1*CF1) )/Rl 

)*R3 

3*ST-SP*T4) )/(RA*Xl*CFt J) ) 

P*T3*R3)/(RA*SP*R3*CF('-)  ) 

l*Gl*X4*G2*G?*X5*G3*G3*X6*(G4*G4BG5*t^5Wx7*  (G6*G6« 

1*A1*X4*a2*A2*X5*A3*a3*X6*(A4*A4>A5*A5WX7« 

) 
l*Gl*X4*A2*G2*X5*A3*G3*X6*<A4*G49A5*r!5WX7* 

) 

l*Ul*Y4*U2*U2*Y5*U3*U3*Y6*(U4*U4eU5*U5WY7* 

) 

k;)-g(N1,  j,K) 

K)-A(NI, J,K) 
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UlsU<NI, J*1»K)-U(N1, J,K) 

G2  =  (i(NI,  J#K*l)^G(NI,  J,K) 

A2kA<N1, J*K*1)-A<NI, J,K> 

U2=U(NI, J»K*1)-U(NI, J,K) 

XJ=tRA*R2*CFl)*(R2-Rl)/R2 

YJsRA*Rl*CFl 

X4srt4*(R4-R3)/(RA*R4«CF(J)  ) 

Y4s|R3*R3>/(PA*R3*Cr(J)) 

A11§A11*0,5*(X3*RI*G1*X4*G2*G2) 

A22»A22*0,5*(X3*A1*A1-^X4*A2*A2) 

A12»A12*0,5*(X3*A1*G1*X4*A2*G2) 

A33|A33*0,5*(Y3*U1*U1*Y4*U2*U2) 
32   CONflNUE 
31   CONflNUE 

All§All*HS 

A12|A12*HS 

A22)|A22*HS 

A33«A33*HS 

A21i|A12 

DEL§A11«A22«A21*A21 

Ql»iA22*rLUl-A21*FLU2)yDFI. 

Q2*(«A21*FLL1*A11*FLU2)/PEL 

Q3«PIU3/A33 

ENES»0.5*(Ql*FLUl*O2*rLU2*O3*FLU3)*0.75*GRA*VOL 

RETURN 

END 


C 
C 
C 


SUBM 

COMP 

USES 

BOUN 

COMH 

liXll 

2B2(^ 

3C3(^ 

4E4(S» 

dGB4{ 

6PEFI 

COMH 

DIMfc 

1GE2I 

2RC1I 

3E!>(i 

NlsN 

N2  =  N 


OUTINE  G 
UTES  THE 
STEEPES 
DARY 

ON  G(9,5 
52),X2(5 
2>.B3(52 
2),C4(52 
2),GX1(5 
52),CF(5 
iBETEiGR 
ON  EX(52 
NSION  DI 
52),6E3( 
2),RC2(2 
),E6<3) 
J*l 
J*2 


RAn(RMAX.PR0l#Qli02»03) 

SQUARE  OF  THE  GRADlENf  ON  THE 
T  PESCENT  METHOD  TO  UPDATE  THE 


FPEE  BOONDAPYiAND 
PfiSITtOW  OF  THE 


?i52> 
2)  1X3 
)iR4( 

)  .ni( 

2>.GX 
?),D5 
A,H,H 
,52), 
F(52, 
52), G 
),RC3 
.E7(3 


,A(g,5 
(52), X 
52).B5 
52).D2 
2(52), 
(52), D 
S.VOL, 
EF{52, 
52).GD 
E4(52) 
(2). PC 
).E8(3 


2.52) 
4(52) 
(52)/ 
(52), 
GX3(5 
6(52) 

ener, 

52)iE 
1(52) 
,RAi( 
4(2), 
).GC1 


,U(9'.52,52),R(52;52).FI52,«2).T(52'.52) 

•Y1(52},Y2(52),Y3(52).V4(5?),B1(52), 

B6(52)»B7(52),B8i52)«Ci(52),C2(52), 

D3(52)'D4(52),El|52),Ea(52»,E3(52)# 

2).Gx4(52),G81(52),GP2J52);gB3(52). 

,D7(52',D8(52),RA,R0.RT,Nr.NJ,PETE» 

RELB 

T(52'.52)'.A11#A12;A22.AJ3 

#GD2(52)'.GD3(52);GC4(52),GP1(52), 

2),RA2(2),Vl(3),V2(3),U3(3),V4(3)i 

REl(2)'»RE2(2),C5l3),C6i3),r7(3»,C8(3), 

(52>'.GC2(52).GC3l52)",Ge4(5?) 
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20 


21 


HA=g,5/HS 

HdsO,5/H 

HCs|H*H)/(HS 

HD»H/(2.0*HS 

S1«0,5*HS 

S2«MS 

S3si,o«0,5*h 

S4sl,,0-HS 

BMAX«0,0 

PHOl=0.0 

DO    I    J"2«N1 

RlaSl*(EX(J# 

R2s«A*Rl*CF( 

R3sP{ j,i)*(S 

R4sf ( J,1)*(S 

PAH(HC*R1*R 

iRUiD) 
R1«S1«(EX(J» 
R2«i8A*Ri*CF( 
R4sF( J,2)*<S 
R4sf ( J,2)*(S 
PA2S(HC*R1*R 

1(EX( J,2)-R(w 
R1  =  SA*S3*F"<^ 
PA3«HC*S3*R1 

1T(J,1))/<R1* 
Rl3«A*S3*F"(  J 
PA4iHC*S3*Rl 

ITU|2)  )/(Rl* 
CONTINUE 
CFH0,5*<Cr( 
R1=0.5*(R( J. 
R2b0.5*(R(J, 
RJ=0,5*(EX(w 
R4sU,5*{EX(w 
PdH(RA*Rl*C 
PB3«(RA*R2*C 
PB6iRA*Rl*CF 
P88|RA*R2*CF 
R5»Q,5*{FUi 
R6sg,5*{F( J« 
R7sQ,5*{EF( J 

R8aQ,5*(EF(J 

R9bS1*(R3-R1 

CONTINUE 

PDli(HD*(R5* 

R98S1*(R4-R2 

PU2s(HD*<R6* 


*HS) 
) 


l)-R(J,l))*R(Jil) 
J) 

i-i,n)-si*Er(J#i) 

1-1.0)-S1*ET(J»1) 
2*(l,0*(R3*R3)/(Rl*Rl)*(R4*R4)/(R2*R2)1l)/(FX(  JiD- 

2)-R(J,2> )*R(J#2> 

J) 

1'1.0)-S1*EF( Ji2> 

l-l,0)-Sl*ET(Ji2) 

2*a.0*{R3*R3)/(Rl*Ri)*{R4*R4)/(R2*R2)^>/ 

.2)) 

,1)*CFU) 

*(1,0*(F(J,1)*F(J,1))/(R( J,1>*R( j,l))*JS3*«3*T(J.l)* 

RD) 

,2)*CF(J) 

*(1.0*(F(J,2)*F(J,2))/tR(  J.2)*R(J,2))*IS3««3*HJ.2)* 

Rl)  ) 


J)*CF 

1)*R( 

2)*R( 

il)*E 

,2)*E 

Fl)*( 

Fl)*( 

1 

1 

1)*F( 

.1)*E 
,2)*E 

)*R1 


(J*l)) 

j*i.in 

J*1.2)) 
X( J*l,l)) 
X( J*l,2)) 
R3.R1)/R1 
R4-R2)/R2 


J*l.l)> 
J*1.2)) 
F(J*1.1)) 

F(J*1.2)) 


(Sl.l,0)-S1*R7)*{RA*R9*CF1)  >/R9 

)tR2 
(Sl-1,0)-S1*R8)*(RA*R9*CF1> )/R9 
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7.2 


23 


PD4?- 

Ri«e. 

R2  =  0, 
PCli( 
PC4i( 
R3  =  0, 
R4sO, 
RbeSl 

PE18( 
P62f- 
CONT 

Ri«y 

R2«A 
XKJ 
X2(J 
X3(J 
Rlag 
X4(J 

RlBli 

R2«A 
YKJ 
Y2(J 
Y3(d 

Rl«U 
Y4(J 
Rl«i* 

R3aA 
BKJ 
B2(J 
B3(J 
CONT 
R2  =  U 
B4(J 
Rls(i 

R3sA 
Bi»(J 
B6(J 
B7(J 
R2*U 
Bb(d 

Rl» 
R3s 
Cl( 
C2( 
C3( 


R2  =  U 
C4(J 


{HD*R 

(HD*R 

5*(R( 

5*(EX 

R1*(R 

R1*R1 

5*{T( 

5*(ET 

•(R2 

HD*R5 

(HD*S 

INUE 

(2,J,1 

(2<J.l 

)»PA1* 

)«PA1* 

)8PA1* 

(N1,J, 
)«PA3* 
(2, J, 2 
(2, J, 2 
)aPA2* 
)sPA2* 
)«PA2* 
(NI,J* 
)sPA4* 
(1,J*1 
(l.J*l 
)BPB1* 

)8PB1* 

)«P81* 

INUE 
(Nl, J* 
)  s 

)sPB3* 

)aPB3* 

)»PB3* 

(NI,J* 

)  s 

(1.J.2 

(l.Ji2 

)«PC1* 

)aPCl* 

)sPCl* 

(NI.J. 

)  s 


5*(RA+S3*CF1*R1) )/Rl 
,6*(RA*S3*CF1*R2) )/R2 
J#l )*R( J,2)) 
(J,1)*EX(J,2)) 
:2«R1)  )/(RA*Rl*CF(J) > 
,>/(RA*Rl*Cr( J) ) 
Jil)*T(J.2n 
(J.1)*ET( J,2)) 


i*(p3*(Sl-1.0)-Sl*R4))/(RA*R5*Cr(J>> 
;3*S3*Rl*R3)/(RA*S3*Rl*cr(J) ) 


)«G(1, J»l) 

)-A(l, J,l) 

R1*R1 

R2*R2 

R1*R2 

1)-U(NI-1,J,1) 

R1*R1 

)<iG(l#  J.2) 

)»A(1, J,2) 

R1*R1 

R2*R2 

R1*R2 

2)-U(Nl-l.J,2) 

R1*R1 

,l)-f;(l,  J#l) 

il)*A(l, J,l) 

R1*R1 

R3*R3 

R1*R3 

1»1)-U(N1, J,l> 

PB6*R2*R2 
•2)*G(1,J.2) 

.2)-A(l, J.2) 

R1*R1 
R3*R3 
R1»R3 
Ii?)-U(NI, J,2> 

PB8*R2*R2 
)"G(1, Jil) 
)»A{1, J,l) 
R1*R1 
R3*R3 
R1*R3 
2)-U(Nl. J,l) 

PC4*R2*R2 
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24 


25 


Rlsli(2,  J*l,l)-G(l 

«J. 

.1) 

R2sU(2,J,l)«R(l, J 

♦1. 

.1) 

R3s4(2, J*l,l)-A(l 

,J. 

1) 

R4sA(2,J.l)-A(l,J 

♦  1. 

1) 

CON? INUE 

Dl(g>sPDl*(Rl*Rl- 

R?* 

'R2) 

D2( J)sPDl*(R3*R3- 

R4< 

'R4) 

D3{ J)aPDl*(Rl*R3- 

R2< 

'R4) 

Rl*g(Nl,J*l,l)-U<NI. 

'1,J,1) 

R2«U(NI, J»1)-U(NI 

•1. 

J*l,l) 

D4( J)«PD3*(Rl*Rl» 

R2* 

'R2> 

R1=U(2, J*l,2)-G(l 

#J. 

2) 

R2=a(2.J,2)»G(l,J 

♦1, 

2) 

R3sA(2,J*l,2)-A(l 

#  J. 

2) 

R4cA(2,J,2)"A(l, J 

♦1. 

2) 

Di>(  J)sPD2*(Rl*Rl' 

R2* 

rR2) 

D6( J)3PD2*(R3*R3- 

R4< 

'R4) 

D7(J)»PD2*(Rl*R3- 

R2< 

rR4) 

Rl=U(NI,J*l,?)-U(NI- 

•l#Ji2) 

R2sU<NI,  J»2)'.U(NI 

•1. 

J*1.2) 

D8( J)ePD4*(Rl*Rl- 

R2* 

'R2) 

CONTINUE 

R1«Q<2,J,2>-G(1.J 

.1) 

R2  =  li(2,J,l)«G(l,  J 

,2) 

R3  =  A(2,  Ji2)-«A(1,  J 

»1) 

R4eA{2,  Jil)"»A(l,  J 

«2) 

El( J)aPEl*<Rl*Rl- 

R2* 

rR2) 

E2( J)sPEl*(R3*R3» 

R4* 

rR4) 

Ei{J)sPEl*(Rl*R3- 

R2* 

'R4) 

RX«U<NI, J#2)-U(NI 

•1. 

Jil) 

R28U<NI,  JiD'-UtNl 

•1. 

J. 2) 

E4( J)sPE2*<Rl*Rl- 

R2* 

■R2) 

CONf INUE 

DO    2    K«2.N1 

Xl{N2)aXl{2) 

XKDsXKNl) 

X2(N2)rx2(2) 

X2(l)«X2(Nl) 

Xi(N2)sX3(2) 

X3(J,)»X3{N1) 

X4(N2)=X4(2> 

X4(l)sX4(Nl) 

YX(N2)«Y1(2) 

YKDaYKNl) 

Y2(N2)sY2{2) 

Y2(l)sY2(Nl) 

Yi(N2)rY3(2)    , 

Y3(1)«Y3(N1) 
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Y4{N2)5Y4(2) 
Y4(l)sY4(Nl) 
BKDsBKNl) 
B2<i)»B2(Nl) 
B3(1,)8B3(N1) 
B4(|,)aB4(Nl) 
B5(l)aB5(Nl) 
B6(i)3B6(Nl) 
B7(X)«B7(N1) 
B8<l)aB8(Nl) 
Cl(N2)»Cl(2) 
CKDbCKNI) 

26  CONTINUE 
C2(N2)eC2(2) 
C2(i)«C2(Nl) 
C3(N2)«C3(2) 
C3(t)aC3(Nl) 
C4(N2)«C4(2) 
C4(t)sC4(Nl) 
Dl(l)8Dl(Nl) 
D2(l)sD2(Nl) 
DJ{l)aD3(Nl) 
D4(l)«D4(Nl) 
D5(1)«D5(N1) 
D6(l)sD6(Nl) 
D7(l)sD7(Nl) 
D8(1)»D8(N1) 
E1(N2)«E1(2) 
EKDsEKNl) 
E2(N2)aE2<2) 
E2(l)sE2(Nl) 
E3(N2)sE3(2) 
E3<J,)3E3(N1) 
E4(N2)sE4(2) 
E4(X)8E4(N1) 

27  CONTINUE 
DO  i    Jel,2 

Rl=Sl*(EX(J,3)-R(J,3))*R(Ji3) 
R2shlA*Ri*CF(  J) 

R4=?(J,3)*(Sl-1.0)-Sl*EF(Ji3) 
R4«T(J,3)*{Sl-l,0)*Sl*ET(Ji3) 
RAl|J)i(HC*Rl*R2*(1.0*(R3*R3)/(Rl*'^l)*{R4*R4W(R?*R2)  )  )/ 

1(EX(J|3)-R(J,3)) 

Rl8KA*S3*F( J,3)*CF( J) 

RA2(  J)sHC*S3*Rl*(l,0*(F(Jj3)*r(j'.3)  )/(R(  Ji3)*R(J,3'l  )♦ 
1(S3*S3*T( J.3)*T( J,3) )/(Rl*Rl) ) 

Rlsa(2i  J#3>.r,(l,  J,S) 

R2  =  A{2,J,3>fiA(l,J,3) 
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Vl( J)8RA1( J)*R1*R1 
V2( J)»RA1( J)*R?*R2 
VJ(J)aRAl( J)*R1*R2 
Rl=U(Nl,J,3)-U(Nl-l,J,3) 
3         V4( J)bRA2( J)*R1*R1 

CFliO,5*(CF(l>*cr(2)) 

R1«0,5*(R(1j3)*R(2,3)) 

R2«0,5*(EX(1,3>*EX(2.3)> 

Rai|(RA*Rl*CFl)*(R2-Rl)/Rl 

R4  =  G<l,2.3).r,(l,l,3) 

R6aA(l,2.3).A(l,l,3) 

Pl=HBl*R4*R4 

P2s«Bl*R6*R6 

P3=HB1*R4*R6 

Ra4l5RA*Rl*CFl 

Rt>sg(Nl,2,3).U(Nl,l,3> 

P^s  RR4*R5*R5 

RJ«0,5#(F(l,3)*r(2,3)) 

R4sQ,5*(EF(l,3)*EF(2.3)) 

Ri>«Sl*(R2-Rl)*Rl 

RD13(HD*(R3*(S1*1,0)-S1*P4)*(RA*R5*CF1)>/R5 

28  CONTINUE 
R5Ba(2,2.3>«G(l,l,3) 
R6Bii(2,l,3>»r,(l,2,3) 
R/=A(2,2.3)»A(1,1,3) 
R8=A(2.1,3)»A(l,2,3) 
Z1  =  «D1*(R5*R«)-R6*R6) 
Z23«Dl*(R7*R7-R8*R8) 
Z3sSDl*(R5*R7'R6*R8) 
R02S-(HD*(RA*S3*R1*CF1)»R3)/R1 
R5«U<NI,2#3)-U(NI-1,1,3) 
R6sU<Nl,l,3)-U(NI-l,2,3) 
Z4e(j|D2*(R5*R5«R6*R6) 

29  CONflNUE 
no  4  Jal,2 

R1«0,5*(R(J,2)*R(J,3)) 
R2eO,5*(EX(j,2)*FXU.3)) 
RC1|J)»(R1*(R2-R1>)/(RA*R1*CF( J)) 
RJ«=(i(l,j,3).iR(l,  J,2) 
R5*A(l,J,3)-A(l,J,2) 
Cb(J)aRcl(J)*R3*R3 

C6( J>aRCl( J)*R5*R5 

C7(J)aRCl(J)*R3*R5 

RC4( J)«(Rl*Ri)/(RA*Rl*CF(J)) 

R4aU(NI,J,3)-U(Nl, J,2) 

C8(J)s  RC4(J)*R4*R4 

R3  =  0.5*(T(  J,2,)*T(J,3)) 

R4=Q,5*(ET(w.2)*ET(J.3)> 
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30   CONflNUE 

Rfcl(v))»(HD*R5*(R3*(Sl*l,0)-Sl*R4)  )  /  (  RA*R5*CF  (  J  )  ) 

R5eG(2,  J,3)-r,(l,  J,2) 

R6  =  G(2,J,2)»r,(l,  J,3) 

R7=A(2i Ji3)-A(l, J,2) 

R8=A(2,J.2)»A(1, J,3) 

E5(  J)8RE1(  J)*(R5*R'5-R6*R6) 

E6( J)sREl( J)*(R7*R7-R8*R8) 

E7( J)8RE1( J)*(R5*R7-R6*R8> 

Rfe2( J)a-(HD*S3*S3*R3*Rl)/(RA*S3*Rl*Cr( J>) 

R!>*U(NI,Ji3)-U(NI-l,  J,2) 

R6=U{NI, J,2)-U(NI-1,J,3) 
4    E8(  J)«RE2(  J)*{R5*R'5-R6*R6) 

DO  5  J-2.N1 

SUMI=0.0 

SUM2«0.0 
SUM3sO.O 

SUM4S0.0 

Rl5Sl*<EX( J*l ,K*1)-R(J*1.K*1> )*R(J*1.K*1) 

R2s«A*Ri*CF(J*l) 

RSsfC  J*l,K*l)*(Sl'»l,0)oSl*EF(J*l,K*l) 

R4s?(j*i,K*l)*(Sl''l,0)=Sl*ET(J*l'.K*l) 

QAli{HC*Rl*R2*(1.0*(R3*R3)/(Rl*Rl)*(R4*R4)/(R2*R?)'l  )/ 
l(EX( J*1,K*1)-R(J*1,K*1)) 

Ri=«(2, J*1,K*1).G(1, J*1,K*1) 

R25A{2,J*l,K*l)«A(l,J*l,K*l) 

Vl(4)sOAl*Rl*Rl 

V2<3)=0A1*R2*R2 

Vi(i)aOAl*Rl*R2 

RlsHA*S3*R< J*1.K*1)*CF(J*1) 
3i   CONTINUE 

OA2iHC*S3*Rl*(l,0*(F(  J*1.K*1)*F{  J*1.K*1)  )/(R«J*l«Kil)*f:(J*l»K*l))* 
l(S3*S3*T(J*l.K*l)*T(J*i,K*l)>/(Rl*Rl)) 

R1»U(NI,J*1,K*1)-U(N1-1,J*1,K*1) 

V4(;i)aOA2*Rl*Rl 

SUMI»SUM1*Y1( J)*n,5*(Yl(J*l)*Xl(J)*Yi(J«l)*Vi{2))+9,25*(Xl(J*l)* 

1V1(1)*X1(J-1)*V1(3)> 

SUM2sSUM2*Y2( J)*0,5*(Y2(J*1)*X2(J)*Y2(J*1)*V2(2))*9,25*(X2<J*1)* 
IV2(1.)*X2(J-1)*V2(3)) 

SUM43SUM3*Y3( J)*0,5*(Y3( J*1)*X3(J)*Y3( J"1)*V3{2> )*9 , 25* ( X3 » J*l )♦ 
lV4(a.)*X3(J-l)*V3(3)) 

SUM4eSUM4*Y4( J)*n,5*(Y4( J*1)*X4( J)*Y4(J»l)*V4(2))*g,25*(X4( J*l)* 

IV4(1)*X4(J-1)*V4(3)) 
32      CONTINUE 

GXHJ)bV1<2) 

Vl(l.)aVl(2> 

V1(2>3V1(3) 
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33 


34 


GX2(J) 
V2(l)s 
V2(2)c 
r,X3fJ) 
VJdJs 
V3(2)« 
GX4(J) 

V4{2)s 
R1«U,5 
R2«0,5 

c^iso, 

OBli(R 
R3=U(1 

RS*A(1 

P5=UBl 

P6=UB1 

P7=QB1 

QB49;RA 

R4sU(N 

P8a 

CONTIN 

SUMl^S 

SUM2=S 

SUM3=S 

SUM4=S 

GBl(J) 
P1  =  H5 

GB2(J) 

P2sP6 

Gti3(J) 

P3»P7 

GB4IJ) 

P4aH8 

R3sO,5 

R4=0.5 

R5=S1* 

nDH(H 

R5=G(2 

R6  =  li(2 

R7s4(2 

CONTIN 

R8=*(2 

ZSsUDl 

Z6sgDl 

Z7=QD1 

Rt)  =  S3* 

0D2g-< 


=V2{2) 

V2{2) 

V2(3) 

sV3(2) 

V3(2) 

V3(3) 

av4(2) 

V4(2) 

V4(3) 

♦(R( Ji 

*(EX{w 

5*(Cr( 

A*R1*C 

,J*liK 

,J*1|K 
♦R3*R3 
*R5*R5 
♦R3*R5 
*R1*CF 
I,J*1« 

UE 

UMl* 

UM2* 

UM3* 

UM4* 

sP5 

sP6 

sp7 


K*1)*R( J*1.K*1) ) 
,K*1)*EX(J*1.K*1) ) 
J)*CF( J*l) ) 
F1)*(R2-R1)/R1 
♦  D-Gd,  J,K*1) 
+  1)-A(1, J,K*1  ) 


K*1)«U(NI, J,K*1> 
QR4*R4*R4 


i5(  J)*B5<  J"!)  ♦0.50*{H1(  J)*BU  J-D+Rl^P*  ) 

6(  J)*B6(.J«'l)  ♦0,50*(B2(  J)*B2«  J'1J*P2*P^  ) 

i7(  J)*B7(.J-1)  ♦0,50*(B3(  J)*B3{  J-1)*P3*P7) 

8(  J)*B8{  J^l)  ♦0,50*(B4(J)*B4(J'»1)*P4*PP  ) 


B5(J) 

B 
B 
B8 


aPfl 


♦  (F< j,K*i)*r( j*i,K*i) ) 

,*{EF(J,K*i)*Er(J*l,K*l)) 

»{R2-R1)*R1 

HD*(R3«<Sl-.l,0)-Sl*R4)*<RA*R5*CFl)  )/l 
?,  J*liK*lUG(l,  J,K*1) 
2, J,K*1)-G(1, J*1,K+1) 
2, J*liK*l).A(l, J,K*1) 

'JUE 

!, J,K*1).A(1, J*1,K*1) 

♦(R5*R5-R6*R6) 

♦(R7*R7-R8*R8) 

*(R5*R7-R6*R8) 

(R2-R1).*R1 

HD*R3*(RA*S3*Rl*cri ) )/Rl 
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3& 


36 


R&=U<NI, J*l#K*l)-U(Ni-l, J,K*1) 
R6=U(NI,J,K*l)-U(NI-l,J*l,K*l) 
Z8=yD2*(R5*R5-R6*R6) 


CONTINUE 

G01iJ)sZ5 

Zl-i5 

GU21J)«Z6 

Z2  =  Z6 

GU3(J>bZ7 

GD4tJ)«Z8 
Z4»28 


C!)(3)»0C1*R3*R3 
C6(i)sQCl*R5*R5 
C7(J)soci*R3*R«5 
f3C4!|(Rl*Rl)/(RA*Rl*Cr(  J*l)  ) 
R4=U(NI,J*1,K*1)-U(NI, J*1,K) 

SUMI»SUM1* 

SUM2=SUM2* 

SUMi=SUM3* 

SUM4=SUM4* 

CONTINUE 

GCllJ)aC5<2) 

C5<i)sC5(2) 

Cb(2)«C5(3) 

GC2|J>»C6(2) 

C6<l)sC6(2) 

C6(2)«C6(3> 

GC3U)«C7(2) 

C7(1)«C7(2) 

C7(2)sC7(3) 

GC4(J)«C8(2) 

C8(l,)sC8{2) 

C8(2)aC8(3) 


i-0<«;jaC8(  J) 

R3=g,5*(T(J*l,K)*T(j*l.K*l)) 
R4=0,5*(ET(w*l.K)*ET(J*l.K*l) ) 
R5»SI*(R2-R1)*R1 
QE1)!(HD*R5*(R3*(S1-1.0)«'S1*R4)  )  /  (  RA*R5*CF  (  J*i  >) 
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R5  =  G(2,J*1,K*1),G(1,J4.1,K) 
R6slj(2,  J*1,K)-G(1,  j*i,K*i  ) 
R7xA(2,J*l,K*l),A(l, j+i,K) 
R8«A(2, J*1,K)-A(1, J*1,K*1) 
Et>(3)»0El*(R5*R5-R6*R6) 
E6{3)=QE1*(R7*97-R8*R8) 
E7<3)aQEi*(R5*R7,R6*RQ) 

n62?-{HD*(S3*S3*R3*Rl) ) / ( RA*S3*Rl*Cp  (  j*i ) ) 
R!>sU(Nr,  J*liK*l)-U(NI.l,  J*1,K) 
R6=U<N1,J*1,K).U(NI-1,J+1,K*1) 
Ea(3)aQE2*(RS*R5-R6*R6) 
37   CONflNUE 

SUMJ.sSUMl*E5<2)*El(J)*0.'s*(El(J*i)*Ei(J-l)*E9(l)*E5(3)) 

SUM2  =  SUM2*E6(2)*E2(J)*0,5*(E2(j*i)*E2(J-l)*E6(n*E4(3)i 

SUM3  =  SUM3*E7(2)*E3(J)*0.5*<E3(j*i)*E3(J*1>*E7(1)*e;!(3)» 

SgM4«SUM4*Ee(2)*E4(J)*0.5*(E4(j*i)*E4(J-l)*Ee{l)*E8(3)) 

GEl|J>sE5(2) 

Et»(l)aE5(2) 

ES(2>=E5(3) 

GS2JJ>bE6(2) 

E6{1)»E6(2) 

E6(2)«E6(3) 

Gfc3{J)»E7(2) 

E7(i>sE7(2) 

E7(2)sE7(3) 

G64|J)«E8(2) 

E8(l)sE8{2) 

Ed{2)=E8(3) 

R2  =  H*H*{RA*R(J.K>*Cr(J))*4,0*R(j'.K) 

Rl=S2*(EX(J,K)-R(J,K>) 

R3=R2*R( J,K) 

GHAUl5{Ql*Ql*SUMl*2.0*f3l*02*SUM3*Q2*Q2*SUM2)/Rl 

GHAU2»(Q3*Q3*SUM4)/R3fGRA 

PR0J,  =  PR01*GRAD? 

Dirj J,K)sGRAni-GRAn2 

DAFsDlFC J,K)/GRAD2 

CC=ABS(DAF) 

PMAXsAMAXKEMAX.CC) 
5    CONflNUE 

DO  6  Jb2iN1 

X1(J)=Y1(J) 

Yl(J)sGXl(J) 

X2(J)sY2{J) 

Y2(J)sGX2(J) 

X3(J)5Y3(J) 

Y3(J)aGX3( J) 

X4(J)sY4(J) 

Y4(J)8GX4(J) 
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6 
7. 


10 


Bi( J)=B5( J) 
Pb( J)sGBl( J) 
B2{ J)aB6( J) 
86( J)sG62( J) 
B3(J)sB7(J) 
B7(J)sGB3(J) 
B4(J)sB8(J) 
Ba( J)sGB4( J) 
C1(J)=GC1(J) 
C2(J)«GC2{J) 
C3(J)sGC3( J) 
C4( J)sRC4( J) 

Dl(J)sD5(J) 
DiJ(J)sQDl(j) 

D2(J)sD6( J) 

D6(J)a6D2(J) 

D4(J>»D7( J) 

D7( J)»GD3( J) 

D4(J)sD8( J) 

De( J)aGD4( J) 

El(J)sGEl(J) 

E2(J)sGE2(J) 

E3(J)sGE3( J) 

E4( J)sGE4( J> 

CONTINUE 

CONTINUE 

DO  lO  Ks2/Nl 

DO  1.0  J  =  2.N1 

R(J|K)sR{J,K)-RELB*Dir(J,K) 

RtTURN 

END 


SUBHOUTINE    SlJRF 

COMMUTES  THE  INITIAL  POSITION  OF  THE  FREE  BOONDARY 

COMMON  G(9,5?,92)#A(9,52.52),U(9'.52,52)iR(52;52).Fl52,«2)»T(52»52) 
liCH9,52)»C2(9,52),El(9,52),E2(9'.52).CF(52),Bl(52);B2(«2)#ni(52), 
2D2<&2),RA#RC,R1,NI,NJ,PETEiPEFI,BETe,GRAiH,HS,V0L.ENER;RELB 

COMMON  EX(52,52),EF(52.52).ET(52'.52).All.A12:A22.Ai3 

N2sNJ*2 

DO  X    J«1,N2 

DO  I  K«l,N2 

R(JiK>aB2(J) 

CONTINUE 

RETURN 

END 
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EOF 


SUBHOUTINE  SaRF 

COMPUTES  THE  POSITION  OF  THE  OUTFR  BOUNDARY 

COMMON  G{9,52,52),A(9,52,52),U(9'.52.52).R(52;52),Fi52.'=2).T(52.52) 
liCl(9,52)#C2(9,52),El(9,52),E2(9'.52).CF(32),51(52):B2(«2).ni{5Z), 
2D2{S2),RA#RC,RI,NI,NJ,PETEiPC-Fl,nETE.GRA.H.HS,V0L,ENER'.RELB 

COMMON  EX(52,52),EF(52,52)iET(52'.52),A11.A12;a?2.A33 

N2=NJ*2 

DO  1  J«liN2 

DO  I  K«1,N2 

EX( JiK)sRO 

CONTINUE 

RETURN 

END 


C 

C 


SUBROUTINE  ElI 

COMPUTES  INITIAL  POSITION  OF  THE  FREE  BOUNDARY, A^D  POSITION  OF 
FJXbD  OUTER  BOUNDARY  FOR  AXIALLY  SYMMETRIC  CASE 

COMMON  GG(9,52,52),GX(9,52#52),GY(^,52.52),Gl(52,5?).G?(52i52). 
103(52,52)  #  0(9,52)  #A(9,52),B(9,52).C{9,52).CFi  52  )  ,  F*  (52  )  ,  R  (  52  ) , 
2EX(52>,EF(52),RA.R0,Rr,NI,NJ,PETF,PAFI,BETE,SRA,M,6S,VfL,ENER, 
3RELB 

ACOR'1,0 

ALARsl.O 

N2=NJ*2 

DO  t  J«1,N2 

Fls(J-2)*H 

Xl=COS(FI)/ACOR 

X2sSlN(Fl )/ALAR 

X3=Xl*xi*X2*x2 

EX( J)3RO*SORT(1,0/X3) 

R(Jl«RI*SQRT(1.0/X3) 

RETURN 

END 
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This  report  was  prepared  as  an  account  of 
Government  sponsored  work.   Neither  the 
United  States,  nor  the  Coininlsslon,  nor  any 
person  acting  on  behalf  of  the  Commission: 

A.  Makes  any  warranty  or  representation, 
express  or  Implied,  with  respect  to  the 
accuracy,  completeness,  or  usefulness  of 
the  Information  contained  In  this  report, 
or  that  the  use  of  any  Information, 
apparatus,  method,  or  process  disclosed 
In  this  report  may  not  Infringe  privately 
owned  rights;  or 

B.  Assumes  any  liabilities  with  respect  to 
the  use  of,  or  for  damages  resulting  from 
the  use  of  any  Information,  apparatus, 
method,  or  process  disclosed  in  this 
report. 

As  used  In  the  above,  "person  acting  on  behalf 
of  the  Commission"  includes  any  employee  or 
contractor  of  the  Commission,  or  employee  of 
such  contractor,  to  the  extent  that  such  em- 
ployee or  contractor  of  the  Commission,  or 
employee  of  such  contractor  prepares,  dis- 
seminates, or  provides  access  to,  any  infor- 
mation pursuant  to  his  employment  or  contract 
with  the  Commission,  or  his  employment  with 
such  contractor. 
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